Soil-water coupled analyzer and soil-water coupled analysis method

ABSTRACT

The soil-water coupled analysis program of the invention computes a global L matrix and a modified global L matrix regarding a volume change rate of a soil skeleton over time, a global H matrix regarding a water permeability of soil, a global M matrix regarding a mass, and a global K matrix regarding a tangent stiffness of the soil skeleton, based on input settings of soils, such as clay, intermediate soil, and sand, to respective elements of a soil foundation, input settings of a solid soil model, and input settings of analysis conditions (steps S 140  to S 165 ). The soil-water coupled analysis program formulates a global tangent stiffness equation (simultaneous linear equations) using all these computed matrixes and determines an unknown ‘jerk field’ and a ‘pore water pressure field’ under given boundary conditions, for example, a given deformation condition and a given stress rate condition (step S 170 ). This enables highly-accurate dynamic and static analyses in soil foundations of various soils from sand to intermediate soils and clay.

TECHNICAL FIELD

The present invention relates to a soil-water coupled analyzer and acorresponding soil-water coupled analysis method. More specifically theinvention pertains to the soil-water coupled analyzer and thecorresponding soil-water coupled analysis method of performingdeformation analysis of a soil skeleton.

BACKGROUND ART

One proposed soil skeleton analysis program is designed to enablethree-dimensional soil deformation analysis or three-dimensionalconsolidation analysis (see, for example, Non-Patent Reference 1).Another proposed soil skeleton analysis program uses a mechanicalconstitutive model expressing static behaviors of a soil foundation(see, for example, Non-Patent Reference 2). These soil skeleton analysisprograms are expected to accurately estimate the water permeability, theconsolidation, the stress in the soil foundation, the shear, the earthpressure, and the slope stability.

-   Non-Patent Reference 1: brochure of ‘FEM analysis support system of    soil structures AFIMEX-GT’, Fujitsu FIP Corporation, No. 0040708-2-   Non-Patent Reference 2: homepage of ‘GeoFEM,    Multi-Purpose/Multi-Physics Parallel Finite Element    Simulator/Platform for Solid Earth’, Coastal Development Institute    of Technology, searched on May 13, 2005.

DISCLOSURE OF THE INVENTION

The prior art soil skeleton analysis program is exclusively used for onesoil type, for example, clay or sand, and does not allow accurateanalysis with regard to intermediate soils as sand and clay mixturesgenerally found in reclaimed artificial soil foundations. The naturalsoil foundation is not made of sand alone or clay alone of an identicalparticle size but includes various types of soils, such as sand-richclay and clay-mixed sand. The results of analysis according to the soilskeleton analysis program developed exclusively for sand or according tothe soil skeleton analysis program developed exclusively for clay arethus significantly different from the actual behaviors of the soilfoundation. The prior art soil skeleton analysis programs are notapplicable to a layered soil foundation of sand and clay.

Rapid repeated loadings in uncompressed and undrained conditions causeliquefaction of a loose sand foundation, while relatively slow repeatedloadings in a drained condition cause compaction of the sand foundation.The prior art soil skeleton analysis programs are applicable to analyzethe liquefaction but are not applicable to analyze the compaction.Namely these analysis programs are usable only exclusively for analysisof the liquefaction. The prior art soil skeleton analysis program isexecutable only after identification of a target soil foundation as asand foundation or a clay foundation suitable for the program with otherindexes. The prior art analysis program is designed exclusively fordynamic analysis or for static analysis and is not applicable to staticanalysis following dynamic analysis as in the case ofliquefaction-induced settlement.

In the soil-water coupled analyzer and the corresponding soil-watercoupled analysis method, there is a demand of enabling both static anddynamic analyses of a soil skeleton. Accurate deformation analysis ofthe soil skeleton is demanded in soil foundations of variousintermediate soils between sand and clay, as well as in layered soilfoundations of clay, sand, and intermediate soil. Another need isiterative deformation analyses of the soil skeleton according to thesoil conditions and the loading conditions. There is another demand ofenabling static deformation analysis following dynamic deformationanalysis, as well as dynamic deformation analysis following staticdeformation analysis.

At least part of the above and the other related demands is attained bya soil-water coupled analyzer and a corresponding soil-water coupledanalysis method of the invention having the configurations discussedbelow.

According to one aspect, the present invention is directed to asoil-water coupled analyzer including: a data input unit that inputsdata on conditions of a soil skeleton and data on external forces; andan analysis module that performs deformation analysis of the soilskeleton. The deformation analysis formulates simultaneous motionequations of rate type with respect to a pore water pressure and avelocity of the soil skeleton including a jerk term of the soilskeleton, based on the input data on the conditions of the soil skeletonand the input data on the external forces, and integrates the formulatedsimultaneous motion equations as needed under a given geometric boundarycondition with regard to any of a displacement, a displacement rate, andan acceleration, a given mechanical boundary condition with regard toeither of a stress and a stress rate, and given hydraulic boundaryconditions with regard to either of the pore water pressure and a totalwater head and a flow of pore water, so as to obtain time historyresponses of a velocity field and a pore water pressure filed.

The soil-water coupled analyzer according to one aspect of the inventionintegrates the simultaneous motion equations of rate type with respectto the pore water pressure and the velocity of the soil skeletonincluding the jerk term of the soil skeleton as needed under the givengeometric boundary condition with regard to any of the displacement, thedisplacement rate, and the acceleration, the given mechanical boundarycondition with regard to either of the stress and the stress rate, andthe given hydraulic boundary conditions with regard to either of thepore water pressure and the total water head and the flow of pore water,in order to obtain the time history responses of the velocity field andthe pore water pressure filed. In the significant influence of theacceleration, the deformation analysis of the soil skeleton is performedby dynamic computations. In the little influence of the acceleration, onthe other hand, the deformation analysis of the soil skeleton isperformed by static computations. The soil-water coupled analyzer ofthis aspect thus enables both dynamic and static deformation analyses ofthe soil skeleton.

In one preferable application of the soil-water coupled analyzer of theinvention, the analysis module calculates multiple terms based on theinput data on the conditions of the soil skeleton and the input data onthe external forces and formulates the simultaneous motion equations ofrate type using the multiple calculated terms. The multiple termsinclude at least one of a first term regarding a volume change rate ofthe soil skeleton over time, a second term regarding a waterpermeability of soil, a third term regarding a mass, and a fourth termregarding a tangent stiffness of the soil skeleton. In this application,the analysis module may formulate the simultaneous motion equations ofrate type using all of the first term, the second term, the third term,and the fourth term. This arrangement ensures appropriate deformationanalysis of the soil skeleton.

In another preferable application of the soil-water coupled analyzer ofthe invention, the analysis module adopts a numerical analysistechnique, such as a finite element method or a difference method, toformulate the simultaneous motion equations of rate type. Thisarrangement desirably facilitates the computations.

In the soil-water coupled analyzer of the invention, one example of thedata on the conditions of the soil skeleton includes data on propertiesof multiple elements constituting the soil skeleton and data onrelations between adjacent elements. The data on the properties of themultiple elements may include data on properties of soils. The settingsare enabled for the properties of the respective multiple elementsconstituting the soil skeleton as well as for the relation between eachpair of adjacent elements. Such setting enables accurate deformationanalysis of the soil skeleton in the soil foundations including variousintermediate soils between sand and clay and accurate deformationanalysis of the soil skeleton in the layered soil foundations of clay,sand, and intermediate soils.

In the soil-water coupled analyzer of the invention, one example of thedata on the external forces includes data on a load and a displacementapplied to a soil-water coupled system. Another example of the data onthe external forces includes data on a vibration applied to a soil-watercoupled system. Such setting enables deformation analysis of the soilskeleton with regard to the load and the displacement or deformationanalysis of the soil skeleton with regard to the vibration.

In one preferable application of the soil-water coupled analyzer of theinvention, the analysis module performs the deformation analysis pluraltimes with changing the data on the external forces. In thisapplication, it is preferable that the analysis module performs thedeformation analysis plural times with setting a result of a currentcycle of the analysis prior to the change of the data on the externalforces to an initial state of the data on the conditions of the soilskeleton for a next cycle of the analysis. This arrangement enablesiterative deformation analyses of the soil skeleton according to thesoil conditions and the loading conditions. This arrangement alsoenables static deformation analysis following dynamic deformationanalysis, as well as dynamic deformation analysis following staticdeformation analysis.

According to another aspect, the invention is also directed to asoil-water coupled analysis method that performs deformation analysis ofa soil skeleton. The deformation analysis formulates simultaneous motionequations of rate type with respect to a pore water pressure and avelocity of the soil skeleton including a jerk term of the soilskeleton, based on data on conditions of the soil skeleton and data onexternal forces, and integrates the formulated simultaneous motionequations as needed under a given geometric boundary condition withregard to any of a displacement, a displacement rate, and anacceleration, a given mechanical boundary condition with regard toeither of a stress and a stress rate, and given hydraulic boundaryconditions with regard to either of the pore water pressure and a totalwater head and a flow of pore water, so as to obtain time historyresponses of a velocity field and a pore water pressure filed.

The soil-water coupled analysis method according to another aspect ofthe invention integrates the simultaneous motion equations of rate typewith respect to the pore water pressure and the velocity of the soilskeleton including the jerk term of the soil skeleton as needed underthe given geometric boundary condition with regard to any of thedisplacement, the displacement rate, and the acceleration, the givenmechanical boundary condition with regard to either of the stress andthe stress rate, and the given hydraulic boundary conditions with regardto either of the pore water pressure and the total water head and theflow of pore water, in order to obtain the time history responses of thevelocity field and the pore water pressure filed. In the significantinfluence of the acceleration, the deformation analysis of the soilskeleton is performed by dynamic computations. In the little influenceof the acceleration, on the other hand, the deformation analysis of thesoil skeleton is performed by static computations. The soil-watercoupled analysis method of this aspect thus enables both dynamic andstatic deformation analyses of the soil skeleton.

In one preferable application of this aspect of the invention, thesoil-water coupled analysis method calculates multiple terms based onthe data on the conditions of the soil skeleton and the data on theexternal forces and formulates the simultaneous motion equations of ratetype using the multiple calculated terms. The multiple terms include atleast one of a first term regarding a volume change rate of the soilskeleton over time, a second term regarding a water permeability ofsoil, a third term regarding a mass, and a fourth term regarding atangent stiffness of the soil skeleton. In this application, thesoil-water coupled analysis method may formulate the simultaneous motionequations of rate type using all of the first term, the second term, thethird term, and the fourth term. This arrangement ensures appropriatedeformation analysis of the soil skeleton.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 schematically illustrates the configuration of a soil-watercoupled analyzer 20 in one embodiment of the invention;

FIG. 2 is a flowchart showing a series of processing according to asoil-water coupled analysis program 30;

FIG. 3 shows a soil foundation having four elements expressed by a 2×2matrix as an analysis object in a concrete example;

FIG. 4 shows local nodes in each element;

FIG. 5 shows one concrete example of a file for inputting the settingsof analysis conditions;

FIG. 6 shows one concrete example of a file for inputting the settingsof analysis conditions;

FIG. 7 shows one concrete example of a file for inputting the settingsof soils;

FIG. 8 shows one concrete example of a file for inputting the settingsof a solid soil model;

FIG. 9 shows a modeled flow of pore water under a two-dimensional planestrain condition;

FIG. 10 shows a modeled flow of pore water under an axisymmetriccondition;

FIG. 11 shows a modeled flow of pore water under a three-dimensionalcondition;

FIG. 12 shows one concrete example of an output file;

FIG. 13 is a flowchart showing a series of processing upon requirementfor repeated execution of the soil-water coupled analysis program withresults of the computation set to new analysis conditions or initialconditions;

FIG. 14 is a graph showing a computation result under loading of alow-frequency microvibration on a triaxial sand specimen in a drainedcondition on a boundary;

FIG. 15 is a graph showing a computation result under loading of ahigh-frequency microvibration on the triaxial sand specimen in thedrained condition on the boundary;

FIG. 16 shows variations in delayed settlement behavior of a naturallydeposited clay foundation under embankment loading (long-term continuoussettlement behavior) against the improvement in mass permeability(overall water permeability of the soil foundation) by the sand drainmethod;

FIG. 17 is a flowchart showing a modified series of processing accordingto the soil-water coupled analysis program 30 with pore water set to acompressible fluid; and

FIG. 18 shows time-settlement curves after a start of seismic vibrationimmediately below embankment as an example of computation withconsideration of the compressibility of pore water.

BEST MODES OF CARRYING OUT THE INVENTION

One mode of carrying out the invention is described below in detail as apreferred embodiment with reference to the accompanied drawings.

FIG. 1 schematically illustrates the configuration of a soil-watercoupled analyzer 20 in one embodiment of the invention. The soil-watercoupled analyzer 20 of the embodiment is constructed by installing asoil-water coupled analysis program 30 as application software in ageneral-purpose computer 22. The soil-water coupled analysis program 30includes a data input module 32 that inputs data for setting theanalysis conditions, the soil properties, and a solid soil model, acomputation module 34 that analyzes the mechanical behaviors ofsaturated soil from the input data according to an elasto-plasticconstitutive equation, and a result output module 36 that outputs theresult of the analysis.

FIG. 2 is a flowchart showing a series of processing according to thesoil-water coupled analysis program 30. The processing of the soil-watercoupled analysis program 30 is described below in detail with referenceto this flowchart. For the simplicity of explanation, a soil foundationas an analysis object in a concrete example is assumed to have fourelements expressed by a 2×2 matrix. FIG. 3 shows the four elements inthe concrete example. As shown in FIG. 3, element numbers 1, 2, 3, and 4are allocated respectively to a lower left element, a lower rightelement, an upper left element, and an upper right element. Node numbers1, 2, 3, 4, 5, 6, 7, 8, and 9 are allocated respectively to a lower leftnode, a lower center node, a lower right node, a middle left node, amiddle center node, a middle right node, an upper left node, an uppercenter node, and an upper right node. Local node numbers 1, 2, 3, and 4are allocated to respective local nodes in each element counterclockwisefrom a lower left local node as shown in FIG. 4. Internode numbers 1, 2,3, and 4 are allocated to respective internodes in each elementcounterclockwise from a bottom internode as shown in FIG. 4.

The soil-water coupled analysis program 30 first inputs settings ofanalysis conditions (step S100). The settings of analysis conditionsinput here may include selection of a procedure for numerical timeintegration based on the difference method, an update frequency of time(the number of computation steps), and selection of a computationcondition, such as two-dimensional plane strain, axisymmetric, orthree-dimensional. FIGS. 5 and 6 show concrete examples of files forinputting the settings of analysis conditions. The inputs in the file ofFIG. 5 include the setting of a computation condition selected among‘plane strain’, ‘axisymmetric’, and ‘three-dimensional’, the setting ofa deformation theory selected between ‘small deformation theory’ and‘finite deformation theory’, the setting of an effective stress ratewith objectivity in the finite deformation theory selected between‘Jaumann's rate of Cauchy stress’ and ‘Green-Naghdi's rate’, the settingof convergence or non-convergence, the setting of a constitutiveequation selected between ‘original Cam-clay model with super-subloadingyield surfaces and rotational hardening’ and ‘modified Cam-clay modelwith super-subloading yield surfaces and rotational hardening’, thesetting of state quantity output data selected between ‘element averageof Gauss points’ and ‘value of specified Gauss point’, the setting of anoutput data format selected between ‘Text’ and ‘Binary’, and the settingof consideration or non-consideration of self weight. The inputs in thefile of FIG. 6 include the settings of a current program number, thenumber of computation steps, a computation time DT per computation step(sec/step) (hereafter may be expressed as Δt), the number of divisionsof excavated load in excavation, the number of divisions of distributedload in application of distributed load, an input file for soil materialparameters and initial values, an input file for mesh data,specification of an output file for initial computation conditions usedin execution of a next program, specification of an output file forcoordinates, specification of an output file for state classification(unloading, loading (plastic compression/expansion,hardening/softening)), specification of an output file for a reactiveforce, specification of an output file for a pore water pressure,specification of an output file for a mean effective stress p′ and ashear stress q, specification of an output file for a degree R ofoverconsolidation (=1/OCR: OCR represents an overconsolidation ratio),specification of an output file for a degree R* of structure,specification of an output file for a slope of a hardening/softeningthreshold line, specification of an output file for a plastic volumestrain, specification of an output file for initial conditions,specification of an output file for a slope of a plasticcompression/expansion threshold line, specification of an output filefor an evolution degree of anisotropy, specification of an output filefor an elapsed time, and output of file***.dat (preceding output filespecifications).

The soil-water coupled analysis program 30 subsequently inputs thesettings of soils as clay, intermediate soil, and sand with regard tothe respective elements of the soil foundation (step S110). The settingsof soils include material parameters and initial conditions. FIG. 7shows one concrete example of a file for inputting the settings ofsoils. The inputs in the file of FIG. 7 include the number ofelasto-plastic soil types, the number of elastic body types,specification of a plastic measure adopted in an evolution rule ofstructure, an intercept N of an isotropic normal consolidation line ofremolded soil, a critical state constant M, a compression index λ, aswelling index κ, a Poisson's ratio ν, a coefficient of permeability k,a specific gravity Gs of soil particles, an initial lateral pressurecoefficient K0, an initial degree R0 of overconsolidation, an initialdegree R*0 of structure, a normal consolidation index m, a U* shape, astructural degradation index a, a structural degradation index b, astructural degradation index c, a rotational hardening index br, arotational hardening limit surface mb, an initial degree of anisotropy,an initial void ratio, a height of soil foundation/specimen, an initialoverburden pressure of the soil foundation, an initial consolidationpressure of the specimen, a cell pressure (tensile: positive), andsettings of a unit system (time, force, and length).

The soil-water coupled analysis program 30 then inputs the settings of asolid soil model (step S120). The settings of the solid soil modelinclude conditions for elements of the soil model. FIG. 8 shows oneconcrete example of a file for inputting the settings of the solid soilmodel shown in FIGS. 3 and 4. In the file of FIG. 8, 1×6 data on the1^(st) line successively show the number of elements, the number ofnodes, and the number of elements with application of distributed loadrate in the solid soil model, the number of nodes per element, thenumber of computation steps, and the number of elastic elements. In thefile, 9×5 data on the 2^(nd) to the 10^(th) lines show each of nodenumbers and a type of boundary condition in an x-direction, a value ofthe boundary condition in the x-direction, a type of boundary conditionin a y-direction, and a value of the boundary condition in they-direction with regard to the corresponding node. The types of boundarycondition include coordinate control (0), load rate control (1),displacement rate control (−1), and displacement acceleration control(−2). In the file of FIG. 8, 4×5 data on the 11^(th) to the 14^(th)lines show each of element numbers and global node numbers mapped tolocal node numbers with regard to the corresponding element (setting ofa mapping of the local node number to the global node number). In thefile, 9×3 data on the 15^(th) to the 23^(rd) lines show each of nodenumbers and a coordinate in the x-direction and a coordinate in they-direction of the corresponding node. Data on the 24^(th) line show thenumber of nodes requiring output of the reactive force, and data on the25^(th) line show their node numbers requiring output of the reactiveforce. There are two output directions of the reactive force, thex-direction and the y-direction with regard to each node. For example,output of the reactive force in the x-direction with regard to the node7 is expressed as 7×2−1=13, and output of the reactive force in they-direction with regard to the node 7 is expressed as 7×2=14. In thefile, 2×6 data on the 26^(th) and 27^(th) lines show each of elementnumbers with application of distributed load rate and a side ofapplication (internode of application), a value of the distributed loadrate of its left node in the x-direction, a value of the distributedload rate in its left node in the y-direction, a value of thedistributed load rate of its right node in the x-direction, and a valueof the distributed load rate of its right node in the y-direction withregard to the corresponding element. In the file, 4×6 data on the28^(th) to the 31^(st) lines show each of element numbers and an elementnumber adjacent to the internode 1 of the corresponding element, anelement number adjacent to the internode 2 of the corresponding element,an element number adjacent to the internode 3 of the correspondingelement, an element number adjacent to the internode 4 of thecorresponding element, and a soil type of the corresponding element(settings of a mapping of the internode number to the global elementnumber and the soil type). Data on the 32^(nd) line shows the number ofconditions satisfying a fixed length between two nodes. When there isany condition satisfying the fixed length between two nodes, itscondition number and specification of nodes involved in the conditionare shown on a next line inserted immediately below the data on the32^(nd) line. Data on the 33^(rd) line shows the number of conditionssatisfying a fixed angle between three nodes. Data on the 34^(th) lineshow a condition number satisfying the fixed angle between three nodesand specification of nodes involved in the condition. Data on the35^(th) line shows the number of conditions satisfying a fixed directionof a relative velocity of two nodes specified by relative positionvectors of former two nodes and latter two nodes. Data on the 36^(th)line show a condition number satisfying the fixed direction of therelative velocity of two nodes and specification of nodes involved inthe condition. Data on the 37^(th) line shows the number of conditionssatisfying an equal displacement (velocity). When there is any conditionsatisfying the equal displacement (velocity), its condition number andspecification relating to the condition are shown on a next lineinserted immediately below the data on the 37^(th) line. Data on the38^(th) line successively show the number of regular waves input intothe model and the number of irregular waves input into the model. Dataon the 39^(th) line successively show a type of regular wave vibration,the number of nodes under influence of the regular wave vibration, andthe amplitude and the period of the regular wave vibration. Data on the40^(th) line shows specification of nodes under influence of the regularwave vibration. There are two directions of the regular wave vibration,the x-direction and the y-direction with regard to each node. Forexample, the x-direction of the regular wave vibration with regard tothe node 7 is expressed as 7×2−1=13, and the y-direction of the regularwave vibration with regard to the node 7 is expressed as 7×2=14. Whenthere is any irregular wave input into the model, specification of theirregular wave is shown on a next line inserted immediately below thedata on the 40^(th) line.

After the data input, series of computations are performed with theinput data (steps S130 to S210). The procedure of this embodiment adoptsthe Wilson's θ method for time integration, although another method, forexample, the Newmark's β method, is also applicable. The procedure ofthis embodiment adopts the Tamura's model or the Christian's model ofthe finite element method of allocating the pore water pressure to theelement center for spatial discretization, although another method, forexample, the Sandhu's model of the finite element method of allocatingthe pore water pressure to the node or the mesh-free method, is alsoapplicable. The following description regards computations according tothe Wilson's θ method and the Tamura's model of the finite elementmethod. A ‘total position vector’, a ‘total velocity increment vector’,and a ‘total acceleration increment vector’, and various conditions (forexample, effective stress, pore water pressure,structure/overconsolidation, and anisotropy) at a time t=t+θΔt areestimated according to Equations (1) to (5) given below (step S130).These predictive calculations are performed for the purpose ofconvergence in a shorter time period and are not essential. In the caseof omission of such predictive calculations, the values eventuallydetermined at a time t=t or the initial values at an initial time areset to the values at the time t=t+θΔt. Each set of vector quantities ofall nodes as in the left side or in each term of the right side ofEquation (1) is expressed by a column vector having components as shownin Equation (6). Each set of scalar quantities of all elements as in theleft side or in each term of the right side of Equation (4) is expressedby a column vector having components as shown in Equation (7).

$\begin{matrix}{{\left\{ x^{N} \right\} ❘_{t + {{\theta\Delta}\; t}}} = {{\left\{ x^{N} \right\}{_{t}{+ \left\{ {v^{N}\left( {{\theta\Delta}\; t} \right)} \right\}}}_{t}} + {\frac{1}{2}\left\{ {{\overset{.}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{2} \right\}{_{t}{{+ \frac{1}{6\;}}\left\{ {{\overset{¨}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{3} \right\}}}_{t}}}} & (1) \\{{\left\{ {v^{N}\left( {{\theta\Delta}\; t} \right)} \right\}{_{t + {{\theta\Delta}\; t}}{= \left\{ {v^{N}\left( {{\theta\Delta}\; t} \right)} \right\}}}_{t}} + {\left\{ {{\overset{.}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{2} \right\}{_{t}{{+ \frac{1}{2\;}}\left\{ {{\overset{¨}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{3} \right\}}}_{t}}} & (2) \\{{\left\{ {{\overset{.}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{2} \right\} ❘_{t + {{\theta\Delta}\; t}}} = {\left\{ {{\overset{.}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{2} \right\}{_{t}{+ \left\{ {{\overset{¨}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{3} \right\}}}_{t}}} & (3) \\{{\left\{ u \right\} ❘_{t + {{\theta\Delta}\; t}}} = {\left\{ u \right\}{_{t}{+ \left\{ {\overset{.}{u}\left( {{\theta\Delta}\; t} \right)} \right\}}}_{t}}} & (4) \\{{A❘_{t + {{\theta\Delta}\; t}}} = {A{_{t}{+ \overset{.}{A}}}_{t}\left( {{\theta\Delta}\; t} \right)}} & (5)\end{matrix}${x^(N)}|_(t), {x^(N)}|_(t+θΔt): column vector (total position vector)expressing set of position vectors of all nodes at time t=t or at timet=t+θΔt{ν^(N)(θΔt)}|_(t), {ν^(N)(θΔt)}|_(t+θΔt): column vector (total incrementvector concerning velocity) given as production of θΔt and column vector(total velocity vector) {ν^(N)} expressing set of velocity vectors ofall nodes at time t=t or at time t=t+θΔt{{dot over (ν)}^(N)(θΔt)²}|_(t), {{dot over (ν)}^(N)(θΔt)²}|_(t+θΔt):column vector (total increment vector concerning acceleration) given asthe production of (θΔt)² and column vector (total acceleration vector){{dot over (ν)}^(N)} expressing set of acceleration vectors of all nodesat time t=t or at time t=t+θΔt{{umlaut over (ν)}^(N)(θΔt)³}|_(t): column vector (total incrementvector concerning jerk) given as production of (θΔt)³ and column vector(total jerk vector) {{umlaut over (ν)}^(N)} expressing set of jerkvectors of all nodes at time t=t{{dot over (u)}(θΔt)}|_(t): column vector (total increment vectorconcerning pore water pressure) given as production of θΔt and columnvector (total pore water pressure rate vector) {{dot over (u)}}expressing set of pore water pressure rate vectors of all elements attime t=t{u}|_(t), {u}|_(t+θΔt): column vector (total pore water pressure vector){u} expressing set of pore water pore vectors of all elements at timet=t or at time t=t+θΔtA: state quantity, for example, effective stress tensor T′, rotationalhardening variable tensor β (tensor quantity), surface force vector t(vector quantity), pore water pressure u, degree R* of structure, degreeR of overconsolidation (scalar quantity)A|_(t), A|_(t+θΔt): state quantity A at time t=t or at time t=t+θΔt{dot over (A)}|_(t): rate of state quantity A at time t=t

$\begin{matrix}{\left\{ a^{N} \right\} = \begin{Bmatrix}a_{1}^{1} \\a_{2}^{1} \\a_{2}^{1} \\\vdots \\a_{1}^{j} \\a_{2}^{j} \\a_{r}^{j} \\\vdots \\a_{1}^{NP} \\a_{2}^{NP} \\a_{r}^{NP}\end{Bmatrix}} & (6)\end{matrix}$a^(j) _(r): r-direction component of a^(N) at j-th nodej: node number (j=1, 2, . . . , NP, NP: total number of nodesr: number of elements per node (r=2 for 2D plane strain condition andaxisymmetric condition, r=3 for 3D condition)

$\begin{matrix}{\left\{ \alpha \right\} = \begin{Bmatrix}\alpha^{1} \\\alpha^{2} \\\vdots \\\alpha^{i} \\\vdots \\\alpha^{({{NE} - 1})} \\\alpha^{NE}\end{Bmatrix}} & (7)\end{matrix}$α^(i): value of α at i-th elementi: element number (i=1, 2, . . . , NE, NE: total number of elements)

A global L matrix (global volume change rate matrix) and a modifiedglobal L matrix (modified global volume change rate matrix) are computedrespectively according to Equation (8) and Equation (10) given below asmatrices for converting the displacement rate of each node to the volumechange rate of the soil skeleton over time (step S140). Each element ofthe global L matrix is shown by Equation (9). The global L matrix has NErows and (NP×r) columns. According to the mapping of the global nodenumber to the local node number specified in the setting of the solidsoil model, the numeral ‘0’ is set to a column component having an m-thnode number mismatching with the global node number in an element Lmatrix for an element i shown in Equation (9). A Bv matrix as theconstituent of the element L matrix is expressed by Equation (12) forthe two-dimensional plane strain condition, by Equation (13) for theaxisymmetric condition, and by Equation (14) for the three-dimensionalcondition. The modified global L matrix is created in the same manner asthe global L matrix and has NE rows and (NP×r) columns.

$\begin{matrix}{L = {\begin{bmatrix}L^{e\; 1} \\L^{e\; 2} \\\vdots \\L^{ei} \\\vdots \\L^{e{({{NE} - 1})}} \\L^{eNE}\end{bmatrix}\text{:}\mspace{14mu}{global}\mspace{14mu} L\mspace{14mu}{matrix}\mspace{14mu}\left( {{global}\mspace{14mu}{volume}\;{change}\mspace{14mu}{rate}\mspace{14mu}{matrix}} \right)}} & (8)\end{matrix}$L ^(ei) =[L ^(i) ¹ ₁ L ^(i) ¹ ₂ L ^(i) ¹ _(r) . . . L ^(i) ^(m) ₁ L ^(i)^(m) ₂ L ^(i) ^(m) ₃ . . . L ^(i) ^(p) ₁ L ^(i) ^(p) ₂ L ^(i) ^(p)_(r)]=∫_(νe) _(i) [B _(ν) ]dν: element L matrix for element i (elementvolume change rate matrix)  (9)

m: local node number in element (m=1, 2, . . . , p)

p: number of nodes per element (p=4 for 2D plane strain condition andaxisymmetric condition, p=8 for 3D condition)

νe^(i): volume area occupied by element i

$\begin{matrix}{\mspace{79mu}{{L_{\theta\; n} = \begin{bmatrix}{\gamma_{\theta\; n}^{1}L^{e\; 1}} \\{\gamma_{\theta\; n}^{2}L^{e\; 2}} \\\vdots \\{\gamma_{\theta\; n}^{i}L^{ei}} \\\vdots \\{\gamma_{\theta\; n}^{({{NE} - 1})}L^{e{({{NE} - 1})}}} \\{\gamma_{\theta\; n}^{NE}L^{eNE}}\end{bmatrix}}\left( {{n = 1},2,3} \right)\text{:}\mspace{14mu}{\quad{\quad{{modified}\mspace{14mu}{global}\mspace{14mu} L\mspace{14mu}{matrix}\mspace{11mu}\left( {{modified}\mspace{14mu}{global}\mspace{14mu}{volume}\mspace{14mu}{change}\mspace{14mu}{rate}\mspace{14mu}{matrix}} \right)}}}}} & (10) \\{\mspace{79mu}{{\gamma_{\theta\; 1}^{i} = {{{- \frac{1}{2{\theta\Delta}\; t}}\frac{k_{i}}{g}} + \frac{1}{6}}},\mspace{79mu}{\gamma_{\theta\; 2}^{i} = {1 - {\frac{1}{{\theta\Delta}\; r}\frac{k_{i}}{g}}}},\mspace{79mu}{\gamma_{\theta\; 3}^{i} = {\frac{1}{6}\left( {2 - {\frac{3}{{\theta\Delta}\; t}\frac{k_{i}}{g}}} \right)}}}} & (11)\end{matrix}$θ: parameter in Wilson's θ methodk_(i): coefficient of permeability of element ig: gravitational acceleration[B _(ν) ]=[N _(,1) ¹ N _(,2) ¹ . . . N _(,1) ^(m) N _(,2) ^(m) . . . N_(,1) ^(p) N _(,2) ^(p)](p=4)  (12)N^(m): shape function regarding m-th node of element i

$\begin{matrix}{\left\lbrack B_{v} \right\rbrack = {\begin{bmatrix}{N_{.1}^{1} + {\frac{N^{1}}{x_{1}^{G}}\mspace{14mu} N_{.2}^{1}\mspace{14mu}\ldots\mspace{14mu} N_{.1}^{m}} +} \\{{\frac{N^{m}}{x_{1}^{G}}\mspace{14mu} N_{.2}^{m}\mspace{14mu}\ldots\mspace{14mu} N_{.1}^{p}} + {\frac{N^{p}}{x_{1}^{G}}\mspace{14mu} N_{.2}^{p}}}\end{bmatrix}\left( {p = 4} \right)}} & (13)\end{matrix}$x^(G) ₁: distance in radial direction from rotational axis to Gausspoint[B _(ν) ]=[N _(,1) ¹ N _(,2) ² N _(,3) ¹ . . . N _(,1) ^(m) N _(,2) ^(m)N _(,3) ^(m) . . . N _(,1) ^(p) N _(,2) ^(p) N _(,3) ^(p)](p=8)  (14)

The integrals like Equation (9) are performed according to the Gaussnumeral integration method after conversion of a target analysis area ina global coordinate system into a local coordinate system as shown inEquation (15). Equation (15) shows volume integration for thethree-dimensional condition. The integration is performed withconsideration of a unit thickness in a perpendicular direction to thetarget analysis area for the two-dimensional plane strain condition,while being performed with consideration of the rotational symmetryabout the axis for the axisymmetric condition.

$\begin{matrix}{{\int{\int{\int{{f\left( {x,y,z} \right)}{\mathbb{d}x}{\mathbb{d}y}{\mathbb{d}z}}}}} = {\sum\limits_{i = 1}^{NG}{\sum\limits_{j = 1}^{NG}{\sum\limits_{k = 1}^{NG}{w_{i}w_{j}w_{k}{f\begin{pmatrix}{{x\left( {\xi_{i},\eta_{j},\zeta_{k}} \right)},} \\{{y\left( {\xi_{i},\eta_{j},\zeta_{k}} \right)},} \\{z\left( {\xi_{i},\eta_{j},\zeta_{k}} \right)}\end{pmatrix}}{J^{\prime}\left( {\xi_{i},\eta_{j},\zeta_{k}} \right)}}}}}} & (15)\end{matrix}$f=f (x, y, z): integrand in integral area in o-xyz coordinate system(for 3D condition)J′: Jacobian determinant for conversion from o-xyz into ξηζ localcoordinate system having domain of [−1,1]NG: number of Gauss pointsξ_(i)(i=1, . . . , NG), η_(j)(j=1, . . . , NG), ζ_(k)(k=1, . . . , NG):coordinates of Gauss points on ξ, η, ζ axesw_(i)(i=1, . . . , NG), w_(j)(j=1, . . . , NG), w_(k)(k=1, . . . , NG):weights of respective Gauss points

A global H matrix (global water permeability matrix) representing waterpermeability of the soil is then computed according to Equation (16)given below (step S145). The hydraulic boundary conditions, such asdrained or undrained condition, are given by varying the components ofan element H matrix shown in Equation (17). For example, the numeral ‘0’is set to a component in the element H matrix corresponding to aninternode in the undrained condition. The global H matrix has NE rowsand NE columns (r: the number of components per node). According to themapping of the internode number of the element to the global elementnumber specified in the setting of the solid soil model, the numeral ‘0’is set to a component having an m-th internode number mismatching withthe global element number in the element H matrix for the element ishown in Equation (17). The last component in the element H matrix ofthe element i in Equation (17) is substituted into an i-th columncomponent of the element i (i-th row) in the global H matrix. Eachcomponent on the right side of Equation (17) is expressed by Equation(18). The meanings of the respective symbols are defined in FIG. 9 forthe two-dimensional plane strain condition, in FIG. 10 for theaxisymmetric condition, and in FIG. 11 for the three-dimensionalcondition.

$\begin{matrix}{\mspace{79mu}{H = {\begin{bmatrix}H^{e\; 1} \\H^{e\; 2} \\\vdots \\H^{ei} \\\vdots \\H^{e{({{NE} - 1})}} \\H^{eNE}\end{bmatrix}\text{:}\mspace{14mu}{global}\mspace{14mu} H\mspace{14mu}{matrix}\mspace{14mu}\left( {{global}\mspace{14mu}{water}\mspace{14mu}{permeability}\mspace{14mu}{matrix}} \right)}}} & (16) \\{H^{ei} = {\left\lbrack {{\alpha_{1}^{ei}\mspace{14mu}\ldots\mspace{14mu}\alpha_{m}^{ei}\mspace{14mu}\ldots\mspace{14mu}\alpha_{s}^{ei}}\mspace{14mu} - {\sum\limits_{m = 1}^{s}\alpha_{m}^{ei}}} \right\rbrack\text{:}\mspace{14mu}{element}\mspace{14mu} H\mspace{14mu}{matrix}\mspace{14mu}\left( {{element}\mspace{14mu}{water}\mspace{14mu}{permeability}\mspace{14mu}{matrix}} \right)}} & (17)\end{matrix}$m internode number in element (m=1, 2, . . . , s)s: number of adjacent elements per element (s=4 for 2D plane straincondition and axisymmetric condition, s=6 for 3D condition)

$\begin{matrix}{\alpha_{m}^{ei} = {\frac{\frac{k_{i}{l_{im} \cdot n^{m}}}{l_{im} \cdot l_{im}} \cdot \frac{k_{m}{l_{mm} \cdot n^{m}}}{l_{mm} \cdot l_{mm}}}{\frac{k_{i}{l_{im} \cdot n^{m}}}{l_{im} \cdot l_{im}} + \frac{k_{m}{l_{mm} \cdot n^{m}}}{l_{mm} \cdot l_{mm}}}S^{m}\mspace{14mu}\left( {{no}\mspace{14mu}{summationn}} \right)}} & (18)\end{matrix}$k_(i): coefficient of permeability of element ik_(m): coefficient of permeability of element mn^(m): <2D plane strain condition and axisymmetric condition>

-   -   outward (seen from element i) unit normal vector perpendicular        to boundary side between element i and element m

<3D condition>

-   -   outward (seen from element i) unit normal vector perpendicular        to two diagonal lines on boundary surface between element i and        element m        S^(m): <2D plane strain condition>    -   distance between two nodes on boundary surface between element I        and element m

<axisymmetric condition>

-   -   (distance between two nodes)×(distance r′ from axis of        rotational symmetry to center of gravity of side) on boundary        surface between element i and element m

<3D condision>

-   -   projected area of boundary surface between element i and element        m to plane that passes through center of gravity of four nodes        defining boundary surface and has normal vector n^(m)    -   l_(im): relative position vector from center of gravity of        element i to center of gravity of element m on boundary surface        l_(mm): relative position vector from center of gravity of        element m to on boundary surface to center of gravity of element        m

A total load increment vector and a total flow increment vector at thetime t=t+θΔt are calculated respectively according to Equation (19) andEquation (34) given below (step S150). A first term on the right side ofEquation (19) is related to a nominal stress rate vector and isexpressed by Equation (20). A second term on the right side of Equation(19) is related to the Green-Naghdi's rate and is expressed by Equation(21). An N matrix on the right side of Equation (20) and a B matrix onthe right side of Equation (21) are expressed by Equations (22) and (23)for the two-dimensional plane strain condition, by Equations (24) and(25) for the axisymmetric condition, and by Equations (26) and (27) forthe three-dimensional condition. The term of the nominal stress ratevector in Equation (20) is written as Equation (28) and is given as thesum of a traction force rate vector externally given through a boundaryas a first term on the right side and an increment of a surface forcevector accompanied with a geometry variation of the boundary as a secondterm on the right side. The traction force rate vector of the first termon the right side at the time t=t+θΔt is given by conversion of atraction force rate vector set as a boundary condition at the timet=t+θΔt according to Equation (29) in conformity with the Wilson's 0method. An n matrix included in the second term on the right side andrelated to an outward unit normal vector on a mechanical boundary ofstress or stress rate is given by Equation (30) for the two-dimensionalplane strain condition and the axisymmetric condition and by Equation(31) for the three-dimensional condition. A column vector including aneffective stress component in a first term on the right side of Equation(21) is given by Equation (32) for the two-dimensional plane straincondition and the axisymmetric condition and by Equation (33) for thethree-dimensional condition.

$\begin{matrix}\begin{matrix}{{\left\{ {\overset{.}{f}{\theta\Delta}\; t} \right\} ❘_{t + {{\theta\Delta}\; t}}} = {{\sum\limits_{j = 1}^{NS}{\,^{\prime}\left\{ {{\overset{.}{f}}^{j}{\theta\Delta}\; t} \right\}}}❘_{t + {{\theta\Delta}\; t}} +}} \\{{\sum\limits_{i = 1}^{NE}{\,^{''}\left\{ {{\overset{.}{f}}^{ei}{\theta\Delta}\; t} \right\}}}❘_{t + {{\theta\Delta}\; t}}{\text{:}\mspace{14mu}{total}\mspace{14mu}{load}}} \\{{increment}\mspace{14mu}{vector}\mspace{14mu}{at}\mspace{14mu}{time}\mspace{14mu} t} \\{= {t + {{\theta\Delta}\; t}}}\end{matrix} & (19)\end{matrix}$

$\sum\limits_{j = 1}^{NS}{\,^{\prime}\text{:}}$operation of adding value of node on j-th side (mechanical boundary ofstress or stress rate) to value of corresponding global node, based onthe mapping of local node number to global node numberNS: number of planes on mechanical boundary (for 3D condition) or numberof sides on mechanical boundary (for 2D plane strain condition andaxisymmetric condition)

${\sum\limits_{i = 1}^{NE}}^{''}\text{:}$operation of adding value of each node in i-th element to value ofcorresponding global node, based on the mapping of local node number toglobal node number{{dot over (f)} ^(j) θΔt}| _(t+θΔt)=∫_(a) [N] ^(T) {{dot over (s)}_(t)}|_(t+θΔt) da×θΔt  (20){{dot over (f)} ^(ei) θΔt}| _(t+θΔt)==∫_(ν) [B] ^(T) {T′ _(Ω)}|_(t+θΔt)dν×θΔt  (21)

$\begin{matrix}{\lbrack N\rbrack = {\begin{bmatrix}N^{1} & 0 & \ldots & N^{m} & 0 & \ldots & N^{p} & 0 \\0 & N^{1} & \ldots & 0 & N^{m} & \ldots & 0 & N^{p}\end{bmatrix}\left( {p = 4} \right)}} & (22) \\{\lbrack B\rbrack = \begin{bmatrix}N_{,1}^{1} & 0 & \ldots & N_{,1}^{m} & 0 & \ldots & N_{,1}^{p} & 0 \\0 & N_{,2}^{1} & \ldots & 0 & N_{,2}^{m} & \ldots & 0 & N_{,2}^{p} \\N_{,2}^{1} & N_{,1}^{1} & \ldots & N_{,2}^{m} & N_{,1}^{m} & \ldots & N_{,2}^{p} & N_{,1}^{p}\end{bmatrix}} & (23) \\{\lbrack N\rbrack = {\begin{bmatrix}N^{1} & 0 & \ldots & N^{m} & 0 & \ldots & N^{p} & 0 \\0 & N^{1} & \ldots & 0 & N^{m} & \ldots & 0 & N^{p}\end{bmatrix}\left( {p = 4} \right)}} & (24) \\{\lbrack B\rbrack = \begin{bmatrix}N_{,1}^{1} & 0 & \ldots & N_{,1}^{m} & 0 & \ldots & N_{,1}^{p} & 0 \\0 & N_{,2}^{1} & \ldots & 0 & N_{,2}^{m} & \ldots & 0 & N_{,2}^{p} \\N_{,2}^{1} & N_{,1}^{1} & \ldots & N_{,2}^{m} & N_{,1}^{m} & \ldots & N_{,2}^{p} & N_{,1}^{p} \\\frac{N^{1}}{x_{1}^{G}} & 0 & \ldots & \frac{N^{m}}{x_{1}^{G}} & 0 & \ldots & \frac{N^{p}}{x_{1}^{G}} & 0\end{bmatrix}} & (25) \\{\lbrack N\rbrack = {\begin{bmatrix}N^{1} & 0 & 0 & \ldots & N^{m} & 0 & 0 & \ldots & N^{p} & 0 & 0 \\0 & N^{1} & 0 & \ldots & 0 & N^{m} & 0 & \ldots & 0 & N^{p} & 0 \\0 & 0 & N^{1} & \ldots & 0 & 0 & N^{m} & \ldots & 0 & 0 & N^{p}\end{bmatrix}\left( {p = 8} \right)}} & (26) \\{\lbrack B\rbrack = \begin{bmatrix}N_{,1}^{1} & 0 & 0 & \ldots & N_{,1}^{m} & 0 & 0 & \ldots & N_{,1}^{p} & 0 & 0 \\0 & N_{,2}^{1} & 0 & \ldots & 0 & N_{,2}^{m} & 0 & \ldots & 0 & N_{,2}^{p} & 0 \\0 & 0 & N_{,3}^{1} & \ldots & 0 & 0 & N_{,3}^{m} & \ldots & 0 & 0 & N_{,3}^{p} \\N_{,2}^{1} & N_{,1}^{1} & 0 & \ldots & N_{,2}^{m} & N_{,1}^{m} & 0 & \ldots & N_{,2}^{p} & N_{,1}^{p} & 0 \\0 & N_{,3}^{1} & N_{,2}^{1} & \ldots & 0 & N_{,3}^{m} & N_{,2}^{m} & \ldots & 0 & N_{,3}^{p} & N_{,2}^{p} \\N_{,3}^{1} & 0 & N_{,1}^{1} & \ldots & N_{,3}^{m} & 0 & N_{,1}^{m} & \ldots & N_{,3}^{p} & 0 & N_{,1}^{p}\end{bmatrix}} & (27)\end{matrix}${{dot over (s)} _(t)}|_(t+θΔt) ×θΔt={{dot over (t)}}| _(t+θΔt) ×θΔt+([B_(ν) ]−[n] _(t+θΔt) [B]){ν^(N) θΔt}| _(t+θΔt) {t}| _(t+θΔt)  (28){{dot over (t)}}| _(t+θΔt) ={{dot over (t)}}| _(t)+θ({{dot over (t)}}|_(t+Δt) −{{dot over (t)}}| _(t))  (29){{dot over (s)}_(t)}|_(t+θΔt): nominal stress rate vector of i-th sidespecified as mechanical boundary at time t=t+θΔt{t}|_(t+θΔt), {{dot over (t)}}|_(t+θΔt): traction force vector acting oni-th side specified as mechanical boundary at time t=t+θΔt and its rate[n]_(t+θΔt): matrix of respective components (direction cosine, n_(r): rdirection component) of outward unit normal vector at i-th side asmechanical boundary[n] _(t+θΔt) =[n ₁ ² n ₂ ² n ₁ n ₂]  (30)[n] _(t+θΔt) =[n ₁ ² n ₂ ² n ₃ ² n ₁ n ₂ n ₂ n ₃ n ₃ n ₁]  (31)

$\begin{matrix}{{\left\{ T_{\Omega}^{\prime} \right\} ❘_{t + {\theta\;\Delta\; t}}} = {\Omega_{c}\begin{Bmatrix}{{- 2}\; T_{12}^{\prime}} \\{2\; T_{12}^{\prime}} \\{T_{11}^{\prime}2\; T_{22}^{\prime}}\end{Bmatrix}}} & (32)\end{matrix}$T′_(ij): i,j components (i=1, 2, 3, j=1, 2, 3) of effective stresstensor T′ at time t=t+θΔtΩ_(c): 2,1 components of material spin tencor Ω={dot over (R)}R^(T) (R:rotation tensor) at time t=t+θΔt axisymmetric condition: value of 4^(th)component=0

$\begin{matrix}{{\left\{ T_{\Omega}^{\prime} \right\} ❘_{t + {\theta\;\Delta\; t}}} = \begin{Bmatrix}\begin{matrix}\begin{matrix}\begin{matrix}\begin{matrix}{2\left( {{\Omega_{2}T_{13}^{\prime}} - {\Omega_{3}T_{12}^{\prime}}} \right)} \\{2\left( {{\Omega_{3}T_{12}^{\prime}} - {\Omega_{1}T_{23}^{\prime}}} \right)}\end{matrix} \\{2\left( {{\Omega_{1}T_{23}^{\prime}} - {\Omega_{2}T_{13}^{\prime}}} \right)}\end{matrix} \\{{{- \Omega_{2}}T_{13}^{\prime}} + {\Omega_{2}T_{23}^{\prime}} + {\Omega_{3}\left( {T_{11}^{\prime} - T_{22}^{\prime}} \right)}}\end{matrix} \\{{\Omega_{1}\left( {T_{22}^{\prime} + T_{34}^{\prime}} \right)} - {\Omega_{2}T_{12}^{\prime}} + {\Omega_{3}T_{13}^{\prime}}}\end{matrix} \\{{\Omega_{1}T_{12}^{\prime}} + {\Omega_{2}\left( {T_{33}^{\prime}T_{11}^{\prime}} \right)} - {\Omega_{3}T_{23}^{\prime}}}\end{Bmatrix}} & (33)\end{matrix}$T′_(ij): i,j components (i=1, 2, 3, j=1, 2, 3) of effective stresstensor T′ at time t=t+θΔtΩ₁, Ω₂, Ω₃: 3,2 components, 1,3 components, and 2,1 components ofmaterial spin tensor Ω={dot over (R)}R^(T) (R: rotation tensor) at timet=t+θΔt

The total flow increment vector at the time t=t+θΔt is expressed byEquation (34) given below and is obtained as the product of a total flowrate vector expressed by Equation (35) and θΔt.{{dot over (f)} _(u)(θΔt)}|_(t+θΔt) ={{dot over (f)} _(u)}|_(t+θΔt)×θΔt: total flow increment vector at time t=t+θΔt  (34)

$\begin{matrix}{{{{\left\{ {\overset{.}{f}}_{u} \right\} ❘_{t + {\theta\;\Delta\; t}}} = {\begin{Bmatrix}\begin{matrix}\begin{matrix}\begin{matrix}\begin{matrix}{{- {\sum\limits_{m = 1}^{s}{\rho_{w}g\;{\alpha_{m}^{e\; 1}\left( {z_{cm}^{1} - z_{c}^{1}} \right)}}}} + q_{1}} \\{{- {\sum\limits_{m = 1}^{s}{\rho_{w}g\;{\alpha_{m}^{e\; 2}\left( {z_{cm}^{2} - z_{c}^{2}} \right)}}}} + q_{2}}\end{matrix} \\\vdots\end{matrix} \\{{- {\sum\limits_{m = 1}^{s}{\rho_{w}g\;{\alpha_{m}^{e\; i}\left( {z_{cm}^{i} - z_{c}^{i}} \right)}}}} + q_{i}}\end{matrix} \\\vdots\end{matrix} \\{{- {\sum\limits_{m = 1}^{s}{\rho_{w}g\;{\alpha_{m}^{e\;{({{NE} - 1})}}\left( {z_{cm}^{({{NE} - 1})} - z_{c}^{({{NE} - 1})}} \right)}}}} + q_{({{NE} - 1})}} \\{{- {\sum\limits_{m = 1}^{s}{\rho_{w}g\;{\alpha_{m}^{eNE}\left( {z_{cm}^{NE} - z_{c}^{NE}} \right)}}}} + q_{NE}}\end{Bmatrix}❘_{t + {\theta\;\Delta\; t}}\text{:}}}{total}\mspace{14mu}{flow}\mspace{14mu}{rate}\mspace{14mu}{vector}\mspace{14mu}{at}\mspace{14mu}{time}\mspace{14mu} t} = {t + {\theta\;\Delta\; t}}} & (35)\end{matrix}$ρ_(w): density of waterz_(c) ^(i): vertical coordinate of center of gravity of element iz_(cm) ^(i): vertical coordinate of center of gravity of element madjacent to element iq_(i): flow of water per unit time externally supplied to element i

A global M matrix (global mass matrix) is computed according to Equation(36) given below (step S160). An element M matrix on the right side ofEquation (37) is expressed by Equation (37). The global M matrix has(NP×r) rows and (NP×r) columns.

$\begin{matrix}{M = {{\sum\limits_{i = 1}^{NE}}^{\prime}{M^{ei}\text{:}{global}\mspace{14mu} M\mspace{14mu}{matrix}\;\left( {{global}\mspace{14mu}{mass}\mspace{14mu}{matrix}} \right)}}} & (36)\end{matrix}$

${\sum\limits_{i = 1}^{NE}}^{\prime}\text{:}$operation of allocating local node number of N matrix [N] at each ofelements from element number i=1 to element number i=NE to global nodenumbers and combining the components from 0, based on the mapping oflocal node number to global node number specified in setting of solidsoil modelM ^(ei)=∫_(ν) ρ[N] ^(T) [N]dν: element M matrix (element mass matrix) ofelement i  (37)ρ: soil density of element i at time

$t = {t + {\theta\;\Delta\;{t\left( {{\rho = {{\frac{e}{1 + e}\rho_{w}} + {\frac{1}{1 + e}\rho_{s}}}},{\rho_{s}\text{:}{density}\mspace{14mu}{of}\mspace{14mu}{soil}\mspace{14mu}{particles}}} \right)}}}$

A global K matrix (a global tangent stiffness matrix of soil skeleton)is computed according to Equation (38) given below (step S165). Theglobal K matrix has (NP×r) rows and (NP×r) columns. An element K matrixas the constituent of the global K matrix is expressed by Equation (39).In Equation (39), a first term on the right side represents a tangentstiffness ensuring mechanical behavior of identified soil material (forexample, sand or clay). A second term on the right side represents acontribution of a time change of soil shape to the tangent stiffness. Athird term on the right side represents a contribution to the tangentstiffness under drastic influence of the acceleration to the soilskeleton, for example, impact load. A fourth term on the right siderepresents a contribution to the tangent stiffness under application ofa body force. Relative differences given to the ‘structural degradationindex’, the ‘normal consolidation index’, and the ‘rotational hardeninglimit constant’ among the material parameters ensure processingcontinuity from sand, intermediate soil, to clay. A D^(ep) matrix(elasto-plastic matrix) included in the first term of Equation (39) isgiven by Equation (40) in the elasto-plastic state (in the loadingstate). Equation (40) is defined by Equations (41) and (42) for thetwo-dimensional plane strain condition, by Equations (43) and (44) forthe axisymmetric condition, and by Equations (45) and (46) for thethree-dimensional condition. The Dep matrix is given by Equation (47) inthe elastic state. The symbols used in Equations (41) through (46) aredefined by Equations (48). The symbols used in Equations (48) are givenby one or combination of Equations (49).

$\begin{matrix}{K = {{\sum\limits_{i = 1}^{NE}}^{\prime}{K^{ei}\text{:}{global}\mspace{14mu} K\mspace{14mu}{matrix}}}} & (38) \\{\left( {{global}\mspace{14mu}{tangent}\mspace{14mu}{stiffness}\mspace{14mu}{matrix}} \right)\mspace{14mu}{of}\mspace{14mu}{soil}\mspace{14mu}{skeleton}} & \; \\{K^{ei} = {{{\int_{v}{{{\lbrack B\rbrack^{T}\left\lbrack D^{ep} \right\rbrack}\lbrack B\rbrack}\ {\mathbb{d}v}}} + {\int_{v}{{{\left\lbrack N^{\prime} \right\rbrack^{T}\left\lbrack T_{1} \right\rbrack}\left\lbrack N^{\prime} \right\rbrack}\ {\mathbb{d}v}}} + {\int{{{\rho_{w}\lbrack N\rbrack}^{T}\lbrack N\rbrack}\left\{ {\overset{.}{v}}^{Ni} \right\}}}}❘_{t + {{\theta\Delta}\; t}}{{\lbrack B\rbrack_{v}{\mathbb{d}v}} - {\int{{\rho_{w}\lbrack N\rbrack}^{T}{\left\{ b \right\}\left\lbrack B_{v} \right\rbrack}{\mathbb{d}v}\text{:}}}}}} & (39)\end{matrix}$element K matrix (element tangent stiffness matix) of element i{{dot over (ν)}^(Ni)}|_(t+θΔt): column vector expressing accelerationvectors of respecfive nodes in element i at time t=t+θΔt{b}: body force vector per unit mass[D ^(ep) ]=[D ₁ ^(ep) ]−[D ₂ ^(ep)]: elasto-plastic matrix  (40)

$\begin{matrix}{\left\lbrack D_{1}^{ep} \right\rbrack = \begin{bmatrix}{a^{\prime} + b^{\prime}} & a^{\prime} & 0 \\a^{\prime} & {a^{\prime} + b^{\prime}} & 0 \\0 & 0 & {b^{\prime}/2}\end{bmatrix}} & (41) \\{\left\lbrack D_{2}^{ep} \right\rbrack = {\frac{1}{e^{\prime}}\begin{bmatrix}\left( {{c^{\prime}{\hat{\eta}}_{11}} - d^{\prime}} \right)^{2} & {\left( {{c^{\prime}{\hat{\eta}}_{11}} - d} \right)\left( {{c^{\prime}{\hat{\eta}}_{22}} - d^{\prime}} \right)} & {c^{\prime}{{\hat{\eta}}_{12}\left( {{c^{\prime}{\hat{\eta}}_{11}} - d^{\prime}} \right)}} \\\; & \left( {{c^{\prime}{\hat{\eta}}_{22}} - d^{\prime}} \right)^{2} & {c^{\prime}{{\hat{\eta}}_{12}\left( {{c^{\prime}{\hat{\eta}}_{22}} - d^{\prime}} \right)}} \\{{sym}.} & \; & \left( {c^{\prime}{\hat{\eta}}_{12}} \right)^{2}\end{bmatrix}}} & (42)\end{matrix}${circumflex over (n)}_(ij): i,j components of tensor {circumflex over(η)}=η−β (i=1, 2, j=1, 2)

$\begin{matrix}{\left\lbrack D_{1}^{ep} \right\rbrack = \begin{bmatrix}{a^{\prime} + b^{\prime}} & a^{\prime} & 0 & a \\a^{\prime} & {a^{\prime} + b^{\prime}} & 0 & a \\0 & 0 & {b^{\prime}/2} & 0 \\a^{\prime} & a^{\prime} & 0 & {a^{\prime} + b^{\prime}}\end{bmatrix}} & (43) \\{\left\lbrack D_{2}^{ep} \right\rbrack = {\frac{1}{e^{\prime}}\begin{bmatrix}\left( {{c^{\prime}{\hat{\eta}}_{11}} - d^{\prime}} \right)^{2} & {\left( {{c^{\prime}{\hat{\eta}}_{11}} - d^{\prime}} \right)\left( {{c^{\prime}{\hat{\eta}}_{22}} - d^{\prime}} \right)} & {c^{\prime}{{\hat{\eta}}_{12}\left( {{c^{\prime}{\hat{\eta}}_{11}} - d^{\prime}} \right)}} & {\left( {{c^{\prime}{\hat{\eta}}_{11}} - d^{\prime}} \right)\left( {{c^{\prime}{\hat{\eta}}_{33}} - d^{\prime}} \right)} \\\; & \left( {{c^{\prime}{\hat{\eta}}_{22}} - d^{\prime}} \right)^{2} & {c^{\prime}{{\hat{\eta}}_{12}\left( {{c^{\prime}{\hat{\eta}}_{22}} - d^{\prime}} \right)}} & {\left( {{c^{\prime}{\hat{\eta}}_{22}} - d^{\prime}} \right)\left( {{c^{\prime}{\hat{\eta}}_{33}} - d^{\prime}} \right)} \\\; & \; & \left( {c^{\prime}{\hat{\eta}}_{12}} \right)^{2} & {c^{\prime}{{\hat{\eta}}_{12}\left( {{c^{\prime}{\hat{\eta}}_{33}} - d^{\prime}} \right)}} \\{{sym}.} & \; & \; & \left( {{c^{\prime}{\hat{\eta}}_{33}} - d^{\prime}} \right)^{2}\end{bmatrix}}} & (44)\end{matrix}${circumflex over (η)}_(ij): i,j components of tensor {circumflex over(η)}=η−β (i=1, 2, 3, j=1, 2, 3)

$\begin{matrix}{\left\lbrack D_{1}^{ep} \right\rbrack = \begin{bmatrix}{a^{\prime} + b^{\prime}} & a^{\prime} & a^{\prime} & \; & \; & \; \\a^{\prime} & {a^{\prime} + b^{\prime}} & a^{\prime} & \; & 0 & \; \\a^{\prime} & a^{\prime} & {a^{\prime} + b^{\prime}} & \; & \; & \; \\\; & \; & \; & {b^{\prime}/2} & \; & \; \\\; & 0 & \; & \; & {b^{\prime}/2} & \; \\\; & \; & \; & \; & \; & {b^{\prime}/2}\end{bmatrix}} & (45) \\{\left\lbrack D_{2}^{ep} \right\rbrack = {\frac{1}{e^{\prime}}\begin{bmatrix}\left( {{c^{\prime}{\hat{\eta}}_{11}} - d^{\prime}} \right)^{2} & {\left( {{c^{\prime}{\hat{\eta}}_{11}} - d^{\prime}} \right)\left( {{c^{\prime}{\hat{\eta}}_{22}} - d^{\prime}} \right)} & {\left( {{c^{\prime}{\hat{\eta}}_{11}} - d^{\prime}} \right)\left( {{c^{\prime}{\hat{\eta}}_{33}} - d^{\prime}} \right)} & {c^{\prime}{{\hat{\eta}}_{12}\left( {{c^{\prime}{\hat{\eta}}_{11}} - d^{\prime}} \right)}} & {c^{\prime}{{\hat{\eta}}_{23}\left( {{c^{\prime}{\hat{\eta}}_{11}} - d^{\prime}} \right)}} & {c^{\prime}{{\hat{\eta}}_{31}\left( {{c^{\prime}{\hat{\eta}}_{11}} - d^{\prime}} \right)}} \\\; & \left( {{c^{\prime}{\hat{\eta}}_{22}} - d^{\prime}} \right)^{2} & {\left( {{c^{\prime}{\hat{\eta}}_{22}} - d} \right)\left( {{c^{\prime}{\hat{\eta}}_{33}} - d^{\prime}} \right)} & {c^{\prime}{{\hat{\eta}}_{12}\left( {{c^{\prime}{\hat{\eta}}_{22}} - d^{\prime}} \right)}} & {c^{\prime}{{\hat{\eta}}_{23}\left( {{c^{\prime}{\hat{\eta}}_{22}} - d^{\prime}} \right)}} & {c^{\prime}{{\hat{\eta}}_{31}\left( {{c^{\prime}{\hat{\eta}}_{22}} - d^{\prime}} \right)}} \\\; & \; & \left( {{c^{\prime}{\hat{\eta}}_{33}} - d^{\prime}} \right)^{2} & {c^{\prime}{{\hat{\eta}}_{12}\left( {{c^{\prime}{\hat{\eta}}_{33}} - d^{\prime}} \right)}} & {c^{\prime}{{\hat{\eta}}_{23}\left( {{c^{\prime}{\hat{\eta}}_{33}} - d^{\prime}} \right)}} & {c^{\prime}{{\hat{\eta}}_{31}\left( {{c^{\prime}{\hat{\eta}}_{33}} - d^{\prime}} \right)}} \\{{{sym}.}\;} & \; & \; & \left( {c^{\prime}{\hat{\eta}}_{12}} \right)^{2} & {\left( {c^{\prime}{\hat{\eta}}_{12}} \right)\left( {c^{\prime}{\hat{\eta}}_{23}} \right)} & {\left( {c^{\prime}{\hat{\eta}}_{12}} \right)\left( {c^{\prime}{\hat{\eta}}_{31}} \right)} \\\; & \; & \; & \; & \left( {c^{\prime}{\hat{\eta}}_{23}} \right)^{2} & {\left( {c^{\prime}{\hat{\eta}}_{23}} \right)\left( {c^{\prime}{\hat{\eta}}_{31}} \right)} \\\; & \; & \; & \; & \; & \left( {c^{\prime}{\hat{\eta}}_{31}} \right)^{2}\end{bmatrix}}} & (46)\end{matrix}${circumflex over (η)}_(ij): i,j components of tensor {circumflex over(η)}=η−β (i=1, 2, 3, j=1, 2, 3)[D₂ ^(ep)]=[0]  (47)

$\begin{matrix}\left. \begin{matrix}{{a^{\prime} = {\overset{\sim}{K} - {\frac{2}{3}\overset{\sim}{G}}}},} \\{{b^{\prime} = {2\;\overset{\sim}{G}}},} \\{{c^{\prime} = {6\;\overset{\sim}{G}}},} \\{{d^{\prime} = {6\;\overset{\sim}{K}\alpha}},} \\{e^{\prime} = {{12\;\eta^{*2}\;\overset{\sim}{G}} + {\overset{\sim}{K}\alpha^{2}} + h}} \\{{\eta\left( {= {q/p^{\prime}}} \right)} = {\sqrt{\frac{3}{2}}{\eta }\left( {= \sqrt{\frac{3}{2}{\eta \cdot \eta}}} \right)\text{:}{effective}\mspace{14mu}{stress}\mspace{14mu}{ratio}}} \\{\beta\text{:}{rotational}\mspace{14mu}{hardening}\mspace{14mu}{variable}\mspace{14mu}{tensor}\;\left( {{tensor}\mspace{14mu}{defining}} \right.} \\{{direction}\mspace{14mu}{and}\mspace{14mu}{magnitute}\mspace{14mu}{of}\mspace{14mu}{anisotropy}\text{:}{The}\mspace{14mu}{greater}} \\{{norm}\;{\beta }{represents}\mspace{14mu}{the}\mspace{14mu}{higher}\mspace{14mu}{degree}\mspace{14mu}{of}\mspace{14mu}{anisotropy}}\end{matrix} \right\} & (48)\end{matrix}$

$\begin{matrix}\left. \begin{matrix}{\overset{\sim}{K} = {\frac{1 + e}{\kappa}{p^{\prime}\left( {= {\frac{J\left( {1 + e_{0}} \right)}{\kappa}p^{\prime}}} \right)}}} \\{\overset{\sim}{G} = {\frac{3\left( {1 - {2\; v}} \right)}{2\left( {1 + v} \right)}\overset{\sim}{K}}} \\{\alpha = {M_{a}^{2} - \eta^{2}}} \\{h = {{Jp}^{\prime}\frac{M^{2} + \eta^{*2}}{MD}\left( {M_{s}^{2} - \eta^{2}} \right)}} \\{D = {\frac{\left( {\lambda - \kappa} \right)}{M\left( {1 + e_{0}} \right)}\text{:}{dilatancy}\mspace{14mu}{factor}}} \\{p^{\prime} = {{- \frac{1}{3}}{trT}^{\prime}\text{:}{mean}\mspace{14mu}{effective}\mspace{14mu}{stress}}} \\{\eta = {{T^{\prime}/p^{\prime}} + I}} \\{\eta^{*} = {\sqrt{\frac{3}{2}}{\hat{\eta}}}}\end{matrix} \right\} & (49) \\{{M_{a}^{2} = {M^{2} + {\frac{3}{2}{\beta \cdot \beta}\text{:}{slope}\mspace{14mu}{of}\mspace{14mu}{threshold}\mspace{14mu}{line}\mspace{14mu}{between}\mspace{14mu}{plastic}}}}{{{compression}\mspace{14mu}{and}\mspace{14mu}{plastic}\mspace{14mu}{expansion}\mspace{14mu}{in}\mspace{14mu} q} - {p^{\prime}\mspace{14mu}{stress}\mspace{14mu}{space}}}} & (50) \\{{M_{s}^{2} = {M_{a}^{2} + {b_{r}\frac{4\; M\;\eta^{*2}}{M^{2} + \eta^{*2}}\left( {{m_{b}\eta^{*}} - {\sqrt{\frac{3}{2}}{\left( {\eta - \beta} \right) \cdot \beta}}} \right)} - {{{MD}\left( {\frac{U^{*}}{R^{*}} - \frac{U}{R}} \right)}\sqrt{{6\;\eta^{*}} + {\frac{1}{3}\alpha^{2}}}\text{:}}}}\begin{matrix}{{slope}\mspace{14mu}{of}\mspace{14mu}{threshold}\mspace{14mu}{line}\mspace{14mu}{between}\mspace{14mu}{hardening}} \\{{{and}\mspace{14mu}{softening}\mspace{14mu}{in}\mspace{14mu} p^{\prime}} - {q\mspace{14mu}{stress}\mspace{14mu}{space}}}\end{matrix}} & (51)\end{matrix}$

-   e: void ratio at time t=t+θΔt-   e₀: void ratio at initial time t=0-   J: determinant of deformation slope tensor of soil skeleton (giving    infinitesimal volume ratio of soil skeleton at time t=t relative to    initial time t=0), J=1 for small deformation analysis-   λ: compression index-   κ: swelling index-   M: critical state constant-   ν: Poisson's ratio-   I: unit tensor-   q=η×p′: shear stress

Equation (50) and Equation (51) respectively show the slope of athreshold line between plastic compression and plastic expansion and theslope of a threshold line between hardening and softening in the p′−qstress space. A plastic stretching norm as the plasticity measure isadopted for the expressions of Equations (50) and (51) as shown inEquations (52) and (53) in the evolution rule, which specifies a processof lowering the degree R* of structure (0<R*≦1: the value of R* closerto 0 represents the more highly structured) and a process of losing thedegree R of overconsolidation (0<R≦1: the value of R closer to 0represents the higher level of overconsolidation) in the structure ofthe soil skeleton. The plasticity measure may be a shear component ofthe plastic stretching or a plastic power given as the inner product ofthe effective stress tensor and the plastic stretching, according to thetype of soil. In small deformation analysis, the plastic stretching isspecified by a plastic strain rate. ‘U*’ in Equation (52) and ‘U’ inEquation (53) are nonnegative functions and are given, for example, byEquations (54) and (55). Equation (56) is adopted for the evolution ruleof the rotational hardening variable tensor representing the anisotropy.

$\begin{matrix}{{\overset{.}{R}}^{*} = {{{JU}^{*}{D^{p}}} = {{JU}^{*}\sqrt{{6\;\eta^{*2}} + {\frac{1}{3}\alpha^{2}}}\frac{{6\overset{\sim}{G}{\hat{\eta} \cdot D}} - {\overset{\sim}{K}{\alpha({trD})}}}{{12\eta^{*2}\overset{\sim}{G}} + {\overset{\sim}{K}\alpha^{2}} + h}}}} & (52) \\{\overset{.}{R} = {{{JU}{D^{p}}} = {{JU}\sqrt{{6\;\eta^{*2}} + {\frac{1}{3}\alpha^{2}}}\frac{{6\overset{\sim}{G}{\hat{\eta} \cdot D}} - {\overset{\sim}{K}{\alpha({trD})}}}{{12\eta^{*2}\overset{\sim}{G}} + {\overset{\sim}{K}\alpha^{2}} + h}}}} & (53)\end{matrix}${dot over (R)}*, {dot over (R)}: rates of R* and R (material timederivation)D^(p): plastic stretchingD: streching

$\begin{matrix}{U^{*} = {\frac{a}{D}{R^{*b}\left( {1 - R^{*}} \right)}^{c}}} & (54) \\{U = {{- \frac{m}{D}}\ln\; R}} & (55)\end{matrix}$a, b, c: structural degradation indexes (material parameters forcontrolling the potential for degrading the structurem: normal consolidation index (material parameters for controlling thepotential for loss of overconsolidation or the potential for normalconsolidation

$\begin{matrix}\begin{matrix}{{\overset{\circ}{\beta}\left( {= {\overset{.}{\beta} + {\beta\;\Omega} - {\Omega\;\beta}}} \right)} = {J\frac{b_{r}}{D}\sqrt{\frac{2}{3}}{D_{s}^{p}}{\hat{\eta}}\left( {{m_{b}\frac{\hat{\eta}}{\hat{\eta}}} - \beta} \right)}} \\{= {J\frac{b_{r}}{D}2\;\eta^{*2}\frac{{6\overset{\sim}{G}{\hat{\eta} \cdot D}} - {\overset{\sim}{K}{\alpha({trD})}}}{{12\eta^{*2}\overset{\sim}{G}} + {\overset{\sim}{K}\alpha^{2}} + h}{\hat{\eta}}\left( {{m_{b}\frac{\hat{\eta}}{\hat{\eta}}} - \beta} \right)}}\end{matrix} & (56)\end{matrix}$

: Green-Nagdhi's rate of β{dot over (β)}: rate of βb_(r): rotational hardening index (material parameter for controllingthe potential for evolution of anisotropy βm_(b): rotational hardening limit constant (material parameter forrepresenting evolution limit of anisotropy

An N′ matrix and a T₁ matrix included in a second term on the right sideof Equation (39) are given by Equations (57) and (58) for thetwo-dimensional plane strain condition, by Equations (59) and (60) forthe axisymmetric condition, and by Equations (61) and (62) for thethree-dimensional condition.

$\begin{matrix}{\left\lbrack N^{\prime} \right\rbrack = \begin{bmatrix}N_{,1}^{1} & 0 & \ldots & N_{,1}^{m} & 0 & \ldots & N_{,1}^{p} & 0 \\0 & N_{,2}^{1} & \ldots & 0 & N_{,2}^{m} & \ldots & 0 & N_{,2}^{p} \\N_{,2}^{1} & 0 & \ldots & N_{,2}^{m} & 0 & \ldots & N_{,2}^{p} & 0 \\0 & N_{,1}^{1} & \ldots & 0 & N_{,1}^{m} & \ldots & 0 & N_{,1}^{p}\end{bmatrix}} & (57) \\{\left\lbrack T_{1} \right\rbrack = \begin{bmatrix}0 & T_{11} & T_{12} & 0 \\T_{22} & 0 & 0 & {- T_{12}} \\T_{12} & 0 & 0 & {- T_{11}} \\0 & T_{12} & {- T_{22}} & 0\end{bmatrix}} & (58)\end{matrix}$T_(ij): i,j components (i=1, 2, j=1, 2) of total stress tensor T(=T′−uI)at time t=t+θΔt (u: pore water pressure, compression: positive)

$\begin{matrix}{\left\lbrack N^{\prime} \right\rbrack = \begin{bmatrix}N_{,1}^{1} & 0 & \ldots & N_{,1}^{m} & 0 & \ldots & N_{,1}^{p} & 0 \\0 & N_{,2}^{1} & \ldots & 0 & N_{,2}^{m} & \ldots & 0 & N_{,2}^{p} \\N_{,2}^{1} & 0 & \ldots & N_{,2}^{m} & 0 & \ldots & N_{,2}^{p} & 0 \\0 & N_{,1}^{1} & \ldots & 0 & N_{,1}^{m} & \ldots & 0 & N_{,1}^{p} \\\frac{N^{1}}{x_{1}^{G}} & 0 & \ldots & \frac{N^{m}}{x_{1}^{G}} & 0 & \ldots & \frac{N^{p}}{x_{1}^{G}} & 0\end{bmatrix}} & (59) \\{\left\lbrack T_{1} \right\rbrack = \begin{bmatrix}0 & T_{11} & {- T_{12}} & 0 & T_{11} \\T_{22} & 0 & 0 & {- T_{12}} & T_{22} \\T_{12} & 0 & 0 & {- T_{11}} & T_{12} \\0 & T_{12} & {- T_{22}} & 0 & T_{12} \\T_{33} & T_{33} & 0 & 0 & 0\end{bmatrix}} & (60)\end{matrix}$T_(ij): i,j components (i=1, 2, 3, j=1, 2, 3) of total stress tensor Tat time t=t+θΔt

$\begin{matrix}{\left\lbrack N^{\prime} \right\rbrack = \begin{bmatrix}N_{,1}^{1} & 0 & 0 & \ldots & N_{,1}^{m} & 0 & 0 & \ldots & N_{,1}^{p} & 0 & 0 \\0 & N_{,2}^{1} & 0 & \ldots & 0 & N_{,2}^{m} & 0 & \ldots & 0 & N_{,2}^{p} & 0 \\0 & 0 & N_{,3}^{1} & \ldots & 0 & 0 & N_{,3}^{1} & \ldots & 0 & 0 & N_{,3}^{p} \\N_{,2}^{1} & 0 & 0 & \ldots & N_{,2}^{m} & 0 & 0 & \ldots & N_{,2}^{p} & 0 & 0 \\0 & N_{,1}^{1} & 0 & \ldots & 0 & N_{,1}^{m} & 0 & \ldots & 0 & N_{,1}^{p} & 0 \\0 & N_{,3}^{1} & 0 & \ldots & 0 & N_{,3}^{m} & 0 & \ldots & 0 & N_{,3}^{p} & 0 \\0 & 0 & N_{,2}^{1} & \ldots & 0 & 0 & N_{,2}^{m} & \ldots & 0 & 0 & N_{,2}^{p} \\0 & 0 & N_{,1}^{1} & \ldots & 0 & 0 & N_{,2}^{m} & \ldots & 0 & 0 & N_{,1}^{p} \\N_{,3}^{1} & 0 & 0 & \ldots & N_{,3}^{m} & 0 & 0 & \ldots & N_{,3}^{p} & 0 & 0\end{bmatrix}} & (61) \\{\left\lbrack T_{1} \right\rbrack = \begin{bmatrix}0 & T_{11} & T_{11} & {- T_{12}} & 0 & 0 & 0 & 0 & {- T_{13}} \\T_{22} & 0 & T_{22} & 0 & {- T_{12}} & {- T_{23}} & 0 & 0 & 0 \\T_{33} & T_{33} & 0 & 0 & 0 & 0 & {- T_{23}} & {- T_{13}} & 0 \\T_{12} & 0 & T_{12} & 0 & {- T_{11}} & {- T_{13}} & 0 & 0 & 0 \\0 & T_{12} & T_{12} & {- T_{22}} & 0 & 0 & 0 & 0 & {- T_{23}} \\T_{23} & T_{23} & 0 & 0 & 0 & 0 & {- T_{22}} & {- T_{12}} & 0 \\T_{23} & 0 & T_{23} & 0 & {- T_{13}} & {- T_{33}} & 0 & 0 & 0 \\0 & T_{31} & T_{31} & {- T_{23}} & 0 & 0 & 0 & 0 & {- T_{33}} \\T_{31} & T_{31} & 0 & 0 & 0 & 0 & {- T_{12}} & {- T_{11}} & 0\end{bmatrix}} & (62)\end{matrix}$T_(ij): i,j components(i=1, 2, 3, j=1, 2, 3) of total stress tensor T attime t=t+θΔt

After calculation of the global M matrix and the global K matrix in theabove manner, an unknown ‘jerk field’ and a pore water pressure fieldare determined according to a global tangent stiffness equation(simultaneous linear equations) given as Equation (63) below under ageometric boundary condition with regard to the displacement, thedisplacement rate, or the acceleration, a mechanical boundary conditionwith regard to the stress or the stress rate, and hydraulic boundaryconditions with regard to the pore water pressure or the total waterhead and the flow of pore water (step S170). Equation (63) is given bycomputing a time derivative term according to an implicit method (finitedifference method) in simultaneous second-order ordinary differentialequations (Equation (64)), which are eventually obtained by finiteelement discretization of the weak form of the motion equations of ratetype and modeling of the soil-water coupled equation with considerationof the influence of an inertial term. Equation (63) represents thesimultaneous linear equations to be (eventually) solved with applicationof the Wilson's θ method to the deformation filed and the trapezoidalrule to the pore water pressure. The total load increment vector, thetotal flow increment vector, the global M matrix, and the global Kmatrix are updated through iterative computations in a time step betweenthe time t=t and the time t=t+θΔt. Equation (64) applies the materialtime derivation of the soil skeleton to the motion equations andaccordingly has a characteristic ‘jerk’ term or second-order derivativeterm on the time of a velocity vector as an explanatory variable. Thischaracteristic gives a rate type constitutive equation of the soilskeleton and allows consideration of a geometric change of the soilfoundation. This enables large deformation analysis requiring accuratemeasurement of a volume change rate of soil between deformation tofailure or with regard to a behavior after the failure. For introductionof this jerk term, Equations (66) to (68) are introduced with respect todeformation based on the concept of a linear acceleration method as thebasis of the Wilson's θ method. In the case of a response having littleeffect of the acceleration, Equation (64) is equal to simultaneousfirst-order ordinary differential equations obtained by static analysis.This means no requirement for the selective use of dynamic analysis orstatic analysis. Under a constraint condition of, for example, a fixedlength between two nodes or a fixed angle between three nodes, imposedon the internode, the Lagrange's method of undetermined multipliers isadopted to solve Equation (64) with consideration of the constraintcondition. The geometric boundary condition on the displacement, thedisplacement rate, or the acceleration set at the time t=t+θΔt isconverted into a jerk vector at the time t=t+θΔt by Equation (69) in thecase of being given as a difference of position vectors or adisplacement vector, by Equation (70) in the case of being given as avelocity vector, and by Equation (71) in the case of being given as anacceleration vector.

$\begin{matrix}{{{\begin{bmatrix}{{\frac{1}{\left( {\theta\;\Delta\; t} \right)^{2}}M} + {\frac{1}{6}K}} & {{- 2}\; L^{T}} \\{- L_{\theta\; 1}} & {H \times \theta\;\Delta\; t}\end{bmatrix}\begin{Bmatrix}{\left\{ {{\overset{¨}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{3} \right\} ❘_{t + {\theta\;\Delta\; t}}} \\{\left\{ u \right\} ❘_{t + {\theta\;\Delta\; t}}}\end{Bmatrix}} = \begin{Bmatrix}\begin{matrix}\begin{matrix}{\left. {\overset{.}{\left\{ f \right.}\left( {\theta\;\Delta\; t} \right)} \right\}{_{t + {\theta\;\Delta\; t}}{{- {K\left\lbrack \left\{ {v^{N}\left( {\theta\;\Delta\; t} \right)} \right\}  \right.}_{t}} + {\left\{ {{\overset{.}{v}}^{N}\left( {\theta\;\Delta\; t} \right)^{2}} \right\}{_{t} +}}}}} \\{\left. {{\frac{1}{3}\left\{ {{\overset{¨}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{3} \right\}}}_{t} \right\rbrack - {2\;{L^{T}\left\lbrack {\left\{ u \right\}{_{t}{{+ \frac{1}{2}}\left\{ {\overset{.}{u}\left( {\theta\;\Delta\; t} \right)} \right\}}}_{t}} \right\rbrack}}}\end{matrix} \\{{\left. {{\overset{.}{\left\{ f \right.}}_{u}\left( {\theta\;\Delta\; t} \right)} \right\}{_{t + {\theta\;\Delta\; t}}{{+ L}\left\{ {v^{N}\left( {\theta\;\Delta\; t} \right)} \right\}}}_{t}} +}\end{matrix} \\{L_{\theta\; 2}\left\{ {{\overset{.}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{2} \right\}{_{t}{{+ L_{\theta\; 3}}\left\{ {{\overset{¨}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{3} \right\}}}_{t}}\end{Bmatrix}}{\left\{ {{\overset{¨}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{3} \right\} ❘_{t + {\theta\;\Delta\; t}}{\text{:}{total}\mspace{14mu}{increment}\mspace{14mu}{vector}\mspace{14mu}{concerning}}}{{{jerk}\mspace{14mu}{at}\mspace{14mu}{time}\mspace{14mu} t} = {t + {\theta\;\Delta\; t}}}} & (63) \\{{{\begin{bmatrix}M & 0 \\0 & 0\end{bmatrix}\frac{\mathbb{d}^{2}}{\mathbb{d}t^{2}}\begin{Bmatrix}\left\{ v^{N} \right\} \\\left\{ u \right\}\end{Bmatrix}} + {\begin{bmatrix}0 & {- L^{T}} \\L^{\prime} & 0\end{bmatrix}\frac{\mathbb{d}}{\mathbb{d}t}\begin{Bmatrix}\left\{ v^{N} \right\} \\\left\{ u \right\}\end{Bmatrix}} + {\begin{bmatrix}K & 0 \\{- L} & H\end{bmatrix}\begin{Bmatrix}\left\{ v^{N} \right\} \\\left\{ u \right\}\end{Bmatrix}}} = \begin{Bmatrix}\left\{ \overset{.}{f} \right\} \\\left\{ {\overset{.}{f}}_{u} \right\}\end{Bmatrix}} & (64) \\{{L^{\prime} = {{\frac{1}{g}\begin{bmatrix}{k_{1}L^{e\; 1}} \\{k_{2}L^{e\; 2}} \\\vdots \\{k_{n}L^{e\; i}} \\\vdots \\{k_{({{NE} - 1})}L^{e{({{NE} - 1})}}} \\{k_{NE}L^{eNE}}\end{bmatrix}}\text{:}}}{{secondary}\mspace{14mu}{modified}\mspace{14mu}{global}\mspace{14mu} L\mspace{14mu}{{matrix}\;\left( {{secondary}\mspace{14mu}{modified}\mspace{14mu}{global}\text{}{volume}\mspace{14mu}{change}\mspace{14mu}{rate}\mspace{14mu}{matrix}} \right)}}} & (65) \\{{\left\{ {\overset{.}{v}}^{N} \right\}{_{t + {\theta\;\Delta\; t}}{= \left\{ {\overset{.}{v}}^{N} \right\}}}_{t}} + {\frac{1}{2}\left\{ {\overset{¨}{v}}^{N} \right\}{_{t}{{\theta\;\Delta\; t} + {\frac{1}{2}\left\{ {\overset{¨}{v}}^{N} \right\}}}}_{t + {\theta\;\Delta\; t}}\theta\;\Delta\; t}} & (66) \\{{\left\{ v^{N} \right\} ❘_{t + {\theta\;\Delta\; t}}} = {{\left\{ v^{N} \right\}{_{t}{+ \left\{ {\overset{.}{v}}^{N} \right\}}}_{t}\theta\;\Delta\; t} + {\frac{1}{3}\left\{ {\overset{¨}{v}}^{N} \right\}{_{t}{\left( {\theta\;\Delta\; t} \right)^{2} + {\frac{1}{6}\left\{ {\overset{¨}{v}}^{N} \right\}}}}_{t + {\theta\;\Delta\; t}}\left( {\theta\;\Delta\; t} \right)^{2}}}} & (67) \\{{\left\{ x^{N} \right\}{_{t + {\theta\;\Delta\; t}}{= \left\{ x^{N} \right\}}}_{t}} + {\left\{ v^{N} \right\}{_{t}{{\theta\;\Delta\; t} + {\frac{1}{2}\left\{ {\overset{.}{v}}^{N} \right\}}}}_{t}\left( {\theta\;\Delta\; t} \right)^{2}} + {\frac{1}{8}\left\{ {\overset{¨}{v}}^{N} \right\}{_{t}{\left( {\theta\;\Delta\; t} \right)^{3} + {\frac{1}{24}\left\{ {\overset{¨}{v}}^{N} \right\}}}}_{t + {\theta\;\Delta\; t}}\left( {\theta\;\Delta\; t} \right)^{3}}} & (68)\end{matrix}${umlaut over (ν)} ^(Nj) ^(r) |_(t+θΔt)(θΔt)³ ={umlaut over (ν)} ^(Nj)^(r) |_(t)(θΔt)³+24θ⁴(x ^(Nj) ^(r) |_(t+Δt) −x ^(Nj) ^(r)|_(t))−24θ³ν^(Nj) ^(r) |_(t)(θΔt)−12θ² {dot over (ν)} ^(Nj) ^(r)|_(t)(θΔt)²−4θ{umlaut over (ν)} ^(Nj) ^(r) |_(t)(θΔt)³  (69){umlaut over (ν)} ^(Nj) ^(r) |_(t+θΔt)(θΔt)³ ={umlaut over (ν)} ^(Nj)^(r) |_(t)(θΔt)³+6θ³(ν^(Nj) ^(r) |_(t+Δt)(θΔt)−ν^(Nj) ^(r)|_(t))(θΔt)−6θ²{dot over (ν)}^(Nj) ^(r) |_(t)(θΔt)²−3θ{umlaut over(ν)}^(Nj) ^(r) |_(t)(θΔt)³  (70){umlaut over (ν)}^(Nj) ^(r) |_(t+θΔt)(θΔt)³={umlaut over (ν)}^(Nj) ^(r)|_(t)(θΔt)³+2θ² {{dot over (ν)} ^(Nj) ^(r) |_(t+Δt)(θΔt)² −{dot over(ν)} ^(Nj) ^(r) |_(t)(θΔt)²}−2θ{umlaut over (ν)} ^(Nj) ^(r)|_(t)(θΔt)³  (71)x^(Nj) ^(r) |_(t|Δt): r-direction component of position vector of j-thnode at time t=t+Δt set as boundary conditionν^(Nj) ^(r) |_(t+Δt): r-direction component of velocity vector of j-thnode at time t=t+Δt set as boundary condition{dot over (ν)}^(Nj) ^(r) |_(t+Δt): r-direction component of accelerationvector of j-th node at time t=t+Δt set as boundary conditionx^(Nj) ^(r) |_(t): r-direction component of position vector of j-th nodeat time t=tν^(Nj) ^(r) |_(t): r-direction component of velocity vector of j-th nodeat time t=t+Δt{dot over (ν)}^(Nj) ^(r) |_(t): r-direction component of accelerationvector of j-th node at time t=t+Δt{umlaut over (ν)}^(Nj) ^(r) |_(t): r-direction component of jerk vectorof j-th node at time t=t+Δt{umlaut over (ν)}^(Nj) ^(r) |_(t+θΔt): r-direction component ofconverted jerk vector of j-th node at time t=t+θΔt

The soil-water coupled analysis program 30 then calculates anacceleration field, a velocity field, coordinates, and various statequantities at the time t=t+θΔt and identifies a loading state of eachGauss point (step S180). A total increment vector concerningacceleration, a total increment vector concerning velocity, and acoordinate vector at the time t=t+θΔt are calculated from the total jerkincrement vector at the time t=t+θΔt in Equation (63) and the variousgeometric quantities at the time t=t according to Equations (72) to(74), which are led from Equations (66) to (68). Respective statequantities A, for example, effective stress, degree R* of structure,degree R of overconsolidation, and anisotropy, at the time t=t+θΔt arecomputed by Equation (75). For example, an effective stress rate at eachGauss point in a certain element at the time t=t+θΔt is calculatedaccording to Equation (76) from velocity vectors of nodes in the certainelement out of the total increment vector concerning velocity given byEquation (73). Subsequent application of Equation (75) determines theeffective stress. The rates of the degree of structure, the degree ofoverconsolidation, and the anisotropy are determined according toEquations (52), (53), and (56) for the Gauss points in theelasto-plastic state and are set equal to 0 for the Gauss points in theelastic state. In the modified Cam-clay model with super-subloadingyield surfaces and rotational hardening, the ‘plastic volume strain’ atthe time t=t+θΔt is calculated by Equation (80) for the Gauss points inthe elasto-plastic state and is kept to a previous value calculated inthe last elasto-plastic state for the Gauss points in the elastic state.In the elastic state, the degree R of overconsolidation at the timet=t+θΔt is computable according to an equation that is obtained bysolving Equation (80) for the degree R of overconsolidation. The porewater pressure at the time t=t+θΔt is directly determined by the globaltangent stiffness equation given as Equation (63). The rate of the porewater pressure is determined by inversely applying Equation (75). Theloading state of each Gauss point in the soil element is identified asthe elasto-plastic state or the elastic state, based on the effectivestress and the stretching of the Gauss point in the soil element. Theloading state of a Gauss point is identified as the elasto-plastic stateupon satisfaction of a loading criterion given by Equation (81) as thenumerator of the plastic multiplier. A column vector of the stretchingis obtained by Equation (82). The components of the column vector at thetime t=t+θΔt are expressed by Equation (83) for the two-dimensionalplane strain condition, by Equation (84) for the axisymmetric condition,and by Equation (85) for the three-dimensional condition.

$\begin{matrix}{{\left\{ {{\overset{.}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{2} \right\}{_{t + {\theta\;\Delta\; t}}{= \left\{ {{\overset{.}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{2} \right\}}}_{t}} + {\frac{1}{2}\left\{ {{\overset{¨}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{3} \right\}{_{t}{{+ \frac{1}{2}}\left\{ {{\overset{¨}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{3} \right\}}}_{t + {\theta\;\Delta\; t}}}} & (72) \\{{\left\{ {v^{N}\left( {\theta\;\Delta\; t} \right)} \right\} ❘_{t + {\theta\;\Delta\; t}}} = {{\left\{ {v^{N}\left( {\theta\;\Delta\; t} \right)} \right\}{_{t}{+ \left\{ {{\overset{.}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{2} \right\}}}_{t}} + {\frac{1}{3}\left\{ {{\overset{¨}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{3} \right\}{_{t}{{+ \frac{1}{6}}\left\{ {{\overset{¨}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{3} \right\}}}_{t + {\theta\;\Delta\; t}}}}} & (73) \\{{\left\{ x^{N} \right\}{_{t + {\theta\;\Delta\; t}}{= \left\{ x^{N} \right\}}}_{t}} + {\left\{ {v^{N}\left( {\theta\;\Delta\; t} \right)}^{2} \right\}{_{t}{{+ \frac{1}{2}}\left\{ {{\overset{.}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{2} \right\}}}_{t}} + {\frac{1}{8}\left\{ {{\overset{¨}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{3} \right\}{_{t}{{+ \frac{1}{24}}\left\{ {{\overset{¨}{v}}^{N}\left( {\theta\;\Delta\; t} \right)}^{3} \right\}}}_{t + {\theta\;\Delta\; t}}}} & (74) \\{{A{_{t + {\theta\;\Delta\; t}}{= A}}_{t}} + {\frac{1}{2}\left( {\overset{.}{A}{_{t + {{\theta\Delta}\; t}}{+ \overset{.}{A}}}_{t}} \right)\left( {\theta\;\Delta\; t} \right)}} & (75)\end{matrix}$A|_(t), A|_(t+θΔt): state quantities A at time t=t and at time t=t+θΔt{dot over (A)}|_(t), {dot over (A)}|_(t+θΔt): rates of state quantitiesA at time t=t and at time t=t+θΔt{{dot over (T)}′}| _(t+θΔt) ×θΔt=[D ^(ep) ][B]{ν ^(Ni) θΔt}| _(t+θΔt)−{T′ _(Ω)}|_(t+θΔt) ×θΔt  (76){ν^(Ni)θΔt}|_(t+Δt): column vector expression of velocity vectors atrespective nodes in element i at time t=t+θΔt{{dot over (T)}′}|_(t+θΔt): column vector expression of components ofeffective stress rate tensor {dot over (T)}′ at respective Gauss pointsin element i at time t=t+θΔt

$\begin{matrix}{{{< {2\; D\mspace{14mu}{plane}\mspace{14mu}{strain}\mspace{14mu}{condition}} > \left\{ {\overset{.}{T}}^{\prime} \right\}}❘_{t + {\theta\;\Delta\; t}}} = {\begin{Bmatrix}\begin{matrix}{\overset{.}{T}}_{11}^{\prime} \\{\overset{.}{T}}_{22}^{\prime}\end{matrix} \\{\overset{.}{T}}_{12}^{\prime}\end{Bmatrix}❘_{t + {\theta\;\Delta\; t}}}} & (77) \\{{{< {{axisymetrix}\mspace{14mu}{condition}} > \left\{ {\overset{.}{T}}^{\prime} \right\}}❘_{t + {\theta\;\Delta\; t}}} = {\begin{Bmatrix}\begin{matrix}{\overset{.}{T}}_{11}^{\prime} \\{\overset{.}{T}}_{22}^{\prime}\end{matrix} \\\begin{matrix}{\overset{.}{T}}_{12}^{\prime} \\{\overset{.}{T}}_{33}^{\prime}\end{matrix}\end{Bmatrix}❘_{t + {\theta\;\Delta\; t}}}} & (78) \\{{{< {3\; D\mspace{14mu}{condition}} > \left\{ {\overset{.}{T}}^{\prime} \right\}}❘_{t + {\theta\;\Delta\; t}}} = {\begin{Bmatrix}\begin{matrix}{\overset{.}{T}}_{11}^{\prime} \\{\overset{.}{T}}_{22}^{\prime}\end{matrix} \\\begin{matrix}{\overset{.}{T}}_{33}^{\prime} \\\begin{matrix}{\overset{.}{T}}_{12}^{\prime} \\{\overset{.}{T}}_{23}^{\prime} \\{\overset{.}{T}}_{31}^{\prime}\end{matrix}\end{matrix}\end{Bmatrix}❘_{t + {\theta\;\Delta\; t}}}} & (79) \\{{{ɛ_{v}^{p}\left( {= {- {\int_{0}^{t}{{JtrD}^{p}\ {\mathbb{d}\tau}}}}} \right)} = {{{MD}\;\ln\frac{p^{\prime}}{p_{0}^{\prime}}} + {{MD}\;\ln\frac{M^{2} + \eta^{*2}}{M^{2}}} + {{MD}\;\ln\; R^{*}} - {{MD}\;\ln\; R\text{:}}}}{{subloading}\mspace{14mu}{yield}\mspace{14mu}{surface}\mspace{14mu}{in}\mspace{14mu}{modified}\mspace{14mu}{Clam}\text{-}{clay}\mspace{14mu}{model}}} & (80)\end{matrix}$ε_(ν) ^(p): ‘plastic volume strain’ (compression: positive)p′₀: p′(>0) value of making the normal yield surface (modified Can-clayyield surface) at the initial time t=0 satisfy η*=0 in the p′−q stressspace

$\begin{matrix}{{{\frac{6\; M\overset{\sim}{G}}{M^{2} + \eta^{*2}}{\left( {\eta - \beta} \right) \cdot D}} - {\overset{\sim}{K}{\alpha({trD})}}} > 0} & (81) \\{{\left\{ D \right\} = {\lbrack B\rbrack\left\{ v^{Ni} \right\}\text{:}}}{{column}\mspace{14mu}{vector}\mspace{14mu}{expression}\mspace{14mu}{of}\mspace{14mu}{components}\mspace{14mu}{of}}\text{}{{stretching}\mspace{14mu}{tensor}\mspace{14mu} D\mspace{14mu}{in}\mspace{14mu}{element}\mspace{14mu} i}} & (82) \\{{{< {2\; D\mspace{14mu}{plane}\mspace{14mu}{strain}\mspace{14mu}{condition}} > \left\{ D \right\}}❘_{t + {\theta\;\Delta\; t}}} = {\begin{Bmatrix}\begin{matrix}D_{11} \\D_{22}\end{matrix} \\{2\; D_{12}}\end{Bmatrix}❘_{t + {\theta\;\Delta\; t}}}} & (83) \\{{{< {{axisymmetric}\mspace{14mu}{condition}} > \left\{ D \right\}}❘_{t + {\theta\;\Delta\; t}}} = {\begin{Bmatrix}\begin{matrix}D_{11} \\D_{22}\end{matrix} \\\begin{matrix}{2\; D_{12}} \\D_{33}\end{matrix}\end{Bmatrix}❘_{t + {\theta\;\Delta\; t}}}} & (84) \\{{{< {3\; D\mspace{14mu}{condition}} > \left\{ D \right\}}❘_{t + {\theta\;\Delta\; t}}} = {\begin{Bmatrix}\begin{matrix}D_{11} \\D_{22}\end{matrix} \\\begin{matrix}D_{33} \\\begin{matrix}{2\; D_{12}} \\{2\; D_{23}} \\{2\; D_{31}}\end{matrix}\end{matrix}\end{Bmatrix}❘_{t + {\theta\;\Delta\; t}}}} & (85)\end{matrix}$

Convergence is then detected according to the effective stresses of therespective Gauss points in the elements (step S190). The procedure ofthis embodiment detects convergence upon satisfaction of Equation (86)for all the Gauss points in all the elements. Here ‘ε’ represents asufficiently small positive value in Equation (86). In the case of nodetection of the convergence, the soil-water coupled analysis program 30returns to step S150 to calculate the total load increment vector andthe total flow increment vector and repeats the subsequent processing ofsteps S150 to S190.

$\begin{matrix}{{\frac{T_{e}^{\prime{(n)}} - T_{e}^{\prime{({n - 1})}}}{T_{e}^{\prime{(n)}}}} \leq ɛ} & (86) \\{{T_{e}^{\prime{(n)}} = {{T^{\prime{(n)}}} = {\sqrt{T^{\prime{(n)}} \cdot T^{\prime{(n)}}}\text{:}}}}{{norm}\mspace{14mu}{of}\mspace{14mu}{effective}\mspace{14mu}{stress}\mspace{14mu}{tensor}\mspace{14mu} T^{\prime{(n)}}\mspace{14mu}{at}\mspace{14mu}{Gauss}\mspace{14mu}{point}}\text{}{{obtained}\mspace{14mu}{by}\mspace{14mu} n\text{-}{times}\mspace{14mu}{iterative}\mspace{14mu}{computations}}\text{}{{{at}\mspace{14mu}{time}\mspace{14mu} t} = {t + {{\theta\Delta}\; t}}}} & (87)\end{matrix}$

In the case of detection of the convergence at step S190, on the otherhand, the soil-water coupled analysis program 30 calculates anacceleration field, a velocity field, coordinates, and various statequantities at a time t=t+Δt and identifies a loading state of each Gausspoint (step S200). The acceleration field, the velocity field, and thecoordinates at the time t=t+Δt are computed from the total incrementvector concerning jerk at the time t=t+θΔt and a total increment vectorconcerning acceleration, a total increment vector concerning velocity,and a coordinate vector at the time t=t according to Equations (88) to(90) given below. Respective state quantities and their rates at thetime t=t+Δt are calculated by Equations (91) and (92) from the statequantities at the time t=t and at the time t=t+θΔt. The loading state ofEach Gauss point at the time t=t+Δt is identified based on therespective state quantities at the time t=t+Δt according to Equations(81) and (82).

$\begin{matrix}{{\left\{ {{\overset{.}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{2} \right\} ❘_{t + {\Delta\; t}}} = {{\left\{ {{\overset{.}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{2} \right\}{_{t}{{+ \frac{1}{\theta}}\left\{ {{\overset{¨}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{3} \right\}}}_{t}} + {\frac{1}{2\;\theta^{2}}\left\lbrack {\left\{ {{\overset{¨}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{3} \right\}{_{t + {{\theta\Delta}\; t}}{- \left\{ {{\overset{¨}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{3} \right\}}}_{t}} \right\rbrack}}} & (88) \\{{\left\{ {v^{N}\left( {{\theta\Delta}\; t} \right)} \right\}{_{t + {\Delta\; t}}{= \left\{ {v^{N}\left( {{\theta\Delta}\; t} \right)} \right\}}}_{t}} + {\frac{1}{\theta}\left\{ {{\overset{.}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{2} \right\}{_{t}{{+ \frac{1}{2\;\theta^{2}}}\left\{ {{\overset{¨}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{3} \right\}}}_{t}} + {\frac{1}{6\;\theta^{3}}\left\lbrack {\left\{ {{\overset{¨}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{3} \right\}{_{t + {{\theta\Delta}\; t}}{- \left\{ {{\overset{¨}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{3} \right\}}}_{t}} \right\rbrack}} & (89) \\{{\left\{ x^{N} \right\} ❘_{t + {\Delta\; t}}} = {{\left\{ x^{N} \right\}{_{t}{{+ \frac{1}{\theta}}\left\{ {v^{N}\left( {{\theta\Delta}\; t} \right)} \right\}}}_{t}} + {\frac{1}{2\;\theta^{2}}\left\{ {{\overset{.}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{2} \right\}{_{t}{{+ \frac{1}{6\;\theta^{3}}}\left\{ {{\overset{¨}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{3} \right\}}}_{t}} + {\frac{1}{24\;\theta^{4}}\left\lbrack {\left\{ {{\overset{¨}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{3} \right\}{_{t + {{\theta\Delta}\; t}}{- \left\{ {{\overset{¨}{v}}^{N}\left( {{\theta\Delta}\; t} \right)}^{3} \right\}}}_{t}} \right\rbrack}}} & (90) \\{\overset{.}{A}{_{t + {\Delta\; t}}{\left( {{\theta\Delta}\; t} \right) = {{2\left( A \right._{t + {{\theta\Delta}\; t}}} - {A\left. _{t} \right)} - \overset{.}{A}}}}_{t}\left( {{\theta\Delta}\; t} \right)} & (91) \\{{A{_{t + {\Delta\; t}}{= A}}_{t}} + {\frac{1}{2\;\theta}\left\{ {\overset{.}{A}{_{t + {\Delta\; t}}{\left( {{\theta\Delta}\; t} \right) + \overset{.}{A}}}_{t}\left( {{\theta\Delta}\; t} \right)} \right\}}} & (92)\end{matrix}$A|_(t+Δt): state quantity A at time t=t+Δt{dot over (A)}|_(t+Δt): rate of state quantity A at time t=t+Δt

It is then determined whether the current computation flow reaches thenumber of computation steps input at step S100 (step S210). Upondissatisfaction of the input number of computation steps, the soil-watercoupled analysis program 30 returns to step S130 to estimate the variousconditions (for example, effective stress, pore water pressure,structure/overconsolidation, and anisotropy) at the time t=t+θΔt. Uponsatisfaction of the input number of computation steps, on the otherhand, the soil-water coupled analysis program 30 is terminated afteroutput of the computation results in the form of a specified output file(step S220).

One concrete example of the output file is shown in FIG. 12. In theoutput file of FIG. 12, data on the 1^(st) line successively show thenumber of elements, the number of nodes, the number of elements withapplication of distributed load rate, the number of nodes withapplication of distributed load rate, the number of nodes per element,the number of nodes×the number of dimensions, the number of nodes perelement×the number of dimensions, the specified number of computationsteps, the total number of implemented computation steps, the number ofdimensions, the number of Gauss points per element, and the number ofelastic elements. Data on the 2^(nd) line successively show an analysiscondition, a deformation theory, an effective stress rate in adoption ofthe finite deformation theory, the setting of convergence ornon-convergence, a constitutive equation, an output method of statequantity data, a format of output data, and the setting of considerationor non-consideration of self weight. Data on the 3^(rd) to the 11^(th)lines successively show each of node numbers, a type of boundarycondition in the x-direction, a value of the boundary condition in thex-direction, a type of boundary condition in the y-direction, and avalue of the boundary condition in the y-direction with regard to thecorresponding node. Data on the 12^(th) line show the number of nodesrequiring output of the reactive force, and data on the 13^(th) lineshow node numbers requiring output of the reactive force as valuesincluding directions at the corresponding nodes. Although being‘omitted’ in the illustrated example of FIG. 12, the output file alsoincludes material parameters at respective Gauss points in each element,values of components at each node, a value of pore water pressure ineach element, values of effective stress components at respective Gausspoints in each element, and a degree of structure. In the output file ofFIG. 12, data on the 1^(st) line after the ‘omission’ show the sequenceof an element number and specified (four) relevant elements of theelement with regard to the elements 1 to 4. Data on the 2^(nd) lineafter the ‘omission’ shows the number of conditions satisfying a fixedlength between two nodes. When the number of conditions is not less than1, its condition number and specification of nodes involved in thecondition are shown on a next line inserted immediately below. Data onthe 3^(rd) line after the ‘omission’ shows the number of conditionssatisfying a fixed angle between three nodes. Data on the 4^(th) lineafter the ‘omission’ show a condition number satisfying the fixed anglebetween three nodes and specification of (three) nodes involved in thecondition. Data on the 5^(th) line after the ‘omission’ shows the numberof conditions satisfying a fixed direction of a relative velocity of twonodes specified by relative position vectors of former two nodes andlatter two nodes. Data on the 6^(th) line after the ‘omission’ show acondition number satisfying the fixed direction of the relative velocityof two nodes and specification of nodes involved in the condition. Dataon the 7^(th) line after the ‘omission’ shows the number of conditionssatisfying an equal displacement (velocity). When there is any conditionsatisfying the equal displacement (velocity), its condition number andspecification relating to the condition are shown on a next lineinserted immediately below. Data on the 8^(th) line after the ‘omission’show the types of soil in the respective elements. Data on the 9^(th)line after the ‘omission’ successively show the number of input regularwaves and the number of input irregular waves. Data on the 10^(th) and11^(th) lines after the ‘omission’ regard the regular waves and show atype of regular wave vibration, the number of nodes under influence ofthe regular wave vibration, the amplitude and the period of the regularwave vibration, and node numbers under influence of the regular wavevibration as values including directions. When there is any inputirregular wave, a type of irregular wave vibration, the number of nodesunder influence of the irregular wave vibration, and node numbers underinfluence of the irregular wave vibration as values including directionsare shown on next two lines. Data on the last line shows a total timesince start of computation.

As described above, the soil-water coupled analysis program 30 of theembodiment inputs multiple different types of soil to be set to therespective elements of the soil foundation. The elasto-plasticdeformation behaviors of the soil skeleton are thus estimated in thesoil-water coupled field under various loading conditions includingdynamic loading conditions and static loading conditions with regard tothe naturally-deposited, highly structured and overconsolidated groundor specimen, as well as the artificially made remolded soil. Thesoil-water coupled analysis program 30 of the embodiment is executed toestimate the deformation behaviors (including bifurcation behavior andpost-peak behavior) and self weight consolidation analysis of thespecimen, the consolidation settlement and lateral movement behaviors(including ‘secondary consolidation’ long-term continuous settlementbehavior) of soil under loading like road embankment, the soildeformation behaviors in excavation, the liquefaction andliquefaction-induced settlement behaviors under repetitive loading likeearthquake, the interaction with an elastic structure (deformationbehaviors of the pile draft foundation), and the deformation behaviorsof reinforced soil. The soil-water coupled analysis program 30 of theembodiment is also executed to study and investigate the waterpumping-induced soil subsidence, the soil bearing capacity, thestability/instability by the groundwater flow, the compacted soilimprovement with the cylindrical cavity expansion by sand piles, thesoil improvement by improvement of mass permeability, and the stabilityof an elasto-plastic embankment structure. In the case of suchestimations and studies, execution of the soil-water coupled analysisprogram plural times may be required with the results of the analysisset for new analysis conditions or initial conditions. In this case, asshown in the flowchart of FIG. 13, after execution of the soil-watercoupled analysis program (step S300), it is determined whether furtherexecution of the program is required (step S310). Upon requirement forfurther execution, the soil-water coupled analysis program is executedagain (step S300) after change of the analysis conditions or the initialconditions (step S320). Upon no requirement for further execution, thesoil-water coupled analysis program is terminated.

Some examples of computation are given with regard to the stiffnessrecovery of the soil skeleton accompanied with compaction andliquefaction-induced settlement. FIG. 14 is a graph showing acomputation result under loading of a low-frequency microvibration on atriaxial sand specimen, and FIG. 15 is a graph showing a computationresult under loading of a high-frequency microvibration on the triaxialsand specimen. As clearly shown by the comparison between these twographs, compaction starts in the course of application of low-frequencymicrovibration load, while liquefaction occurs in the course ofapplication of high-frequency microvibration load and compaction startsafter completion of the load application. Under the condition of a fixedcycle rate, application of the high-frequency microvibration load causesthe greater compaction. The liquefaction-induced settlement is simulatedby first executing the soil-water coupled analysis program under thevibration loading condition in the initial state and subsequentlyexecuting the soil-water coupled analysis program under the vibrationunloading condition in the computation result set to the new initialstate.

As a possible measure against delayed settlement of a naturallydeposited clay foundation, FIG. 16 shows variations in delayedsettlement behavior of the naturally deposited clay foundation underembankment loading (long-term continuous settlement behavior) againstthe improvement in mass permeability (overall water permeability of thesoil foundation) by the sand drain method. The higher mass permeability‘k’ represents the enhanced water permeability. As clearly shown in thisgraph, the unimproved water permeability leads to duration of settlementover 50 years and gives a greater remaining settlement after completionof embankment loading. The improved water permeability, on the otherhand, causes early convergence of settlement and gives a less residualsettlement after completion of embankment loading.

As described in detail above, the soil-water coupled analyzer 20 of theembodiment determines the unknown ‘jerk field’ and the pore waterpressure field according to the global tangent stiffness equation(simultaneous linear equations) under the geometric boundary conditionwith regard to the displacement, the displacement rate, or theacceleration, the mechanical boundary condition with regard to thestress or the stress rate, and the hydraulic boundary conditions withregard to the pore water pressure or the total water head and the flowof pore water. The global tangent stiffness equation is formulated byusing the global L matrix, the modified global L matrix, the global Hmatrix, the global M matrix, and the global K matrix. This allowsconsideration of the geometric change of the soil foundation and enableslarge deformation analysis requiring accurate measurement of the volumechange rate of soil between deformation to failure or with regard to thebehavior after the failure. The soil-water coupled analyzer 20 of theembodiment sets the type of soil for each element of the soil foundationand is thus preferably applicable to computation of mechanical behaviorsin a wide range of soil skeleton structure (structure,overconsolidation, anisotropy) including a continuously varying soilstructure of sand, intermediate soil to clay, unusual soils likevolcanic ash soil, and a layered soil structure of sand, clay, andintermediate soil. In the case of a response having little effect of theacceleration, the global tangent stiffness equation is equal to a globalstiffness equation obtained by static analysis. The mechanical behaviorsof the soil skeleton are thus computable by both dynamic analysis andstatic analysis. The soil-water coupled analyzer 20 of the embodiment isapplicable to study not only the consolidation deformation behaviorsafter earthquake-induced liquefaction but the static deformationanalysis following the dynamic deformation analysis, such as seismicdeformation behaviors during consolidation deformation, and the dynamicdeformation analysis following the static deformation analysis.

The soil-water coupled analyzer 20 of the embodiment sets the type ofsoil for each element of the soil foundation. The possible options ofthe soil type may be limited. The type of soil may be set for each blockof elements, instead of for each element.

The soil-water coupled analyzer 20 of the embodiment accuratelycalculates the global L matrix, the modified global L matrix, the globalH matrix, the global M matrix, and the global K matrix with regard totheir respective element matrixes. The calculation may, however, be notstrict but approximate.

The soil-water coupled analyzer 20 of the embodiment calculates theglobal L matrix, the modified global L matrix, the global H matrix, theglobal M matrix, and the global K matrix and formulates the globaltangent stiffness equation (simultaneous linear equations) by using allthese matrixes. One possible modification may formulate the globaltangent stiffness equation (simultaneous linear equations) by using onlypart of the global L matrix, the modified global L matrix, the global Hmatrix, the global M matrix, and the global K matrix. These matrixes arenot restrictive but any other suitable matrixes may be used for the samepurpose.

The soil-water coupled analyzer 20 of the embodiment adopts the finiteelement method for deformation analysis of the soil skeleton. The finiteelement method is, however, not essential, and any other suitabletechnique may be adopted for deformation analysis of the soil skeleton.

The soil-water coupled analyzer 20 of the embodiment inputs the settingsof analysis conditions as shown in the concrete examples of FIGS. 5 and6. The input analysis conditions include the setting of the computationcondition selected among ‘plane strain’, ‘axisymmetric’, and‘three-dimensional’, the setting of the deformation theory selectedbetween ‘small deformation theory’ and ‘finite deformation theory’, thesetting of the effective stress rate with objectivity in the finitedeformation theory selected between ‘Jaumann's rate of Cauchy stress’and ‘Green-Naghdi's rate’, the setting of convergence ornon-convergence, the setting of the constitutive equation selectedbetween ‘original Cam-clay model with super-subloading yield surfacesand rotational hardening’ and ‘modified Cam-clay model withsuper-subloading yield surfaces and rotational hardening’, the settingof state quantity output data selected between ‘element average of Gausspoints’ and ‘value of specified Gauss point’, the setting of the outputdata format selected between ‘Text’ and ‘Binary’, and the setting ofconsideration or non-consideration of self weight. The input analysisconditions further include the settings of a current program number, thenumber of computation steps, the computation time DT per one computationstep (sec/step), the number of divisions of excavated load inexcavation, the number of divisions of distributed load in applicationof distributed load, the input file for soil material parameters andinitial values, the input file for mesh data, specification of theoutput file for initial computation conditions used in execution of anext program, specification of the output file for coordinates,specification of the output file for state classification (unloading,loading (plastic compression/expansion, hardening/softening)),specification of the output file for the reactive force, specificationof the output file for the pore water pressure, specification of theoutput file for the mean effective stress p′ and the shear stress q,specification of the output file for the degree R of overconsolidation(=1/OCR: OCR represents the overconsolidation ratio), specification ofthe output file for the degree R* of structure, specification of theoutput file for the slope of the hardening/softening threshold line,specification of the output file for the plastic volume strain,specification of the output file for initial conditions, specificationof the output file for the slope of the plastic compression/expansionthreshold line, specification of the output file for the evolutiondegree of anisotropy, specification of the output file for the elapsedtime, and output of file***.dat (preceding output file specifications).The input analysis conditions are, however, not restricted to theseparameters. Only part of these parameters may be input as the analysisconditions or any other suitable parameters may be input as the analysisconditions. These analysis conditions may be input in any suitable inputformat.

The soil-water coupled analyzer 20 of the embodiment inputs the settingsof soils as clay, intermediate soil, and sand with regard to therespective elements of the soil foundation as shown in the concreteexample of FIG. 7. The input settings of soils include the number ofsoil types, the number of elastic body types, specification of theplastic measure adopted in the evolution rule of structure, theintercept N of the isotropic normal consolidation line of remolded soil,the critical state constant M, the compression index λ, the swellingindex κ, the Poisson's ratio ν, the coefficient of permeability k, thespecific gravity Gs of soil particles, the initial lateral pressurecoefficient K0, the initial degree R0 of overconsolidation, the initialdegree R*0 of structure, the normal consolidation index m, the U* shape,the structural degradation index a, the structural degradation index b,the structural degradation index c, the rotational hardening index br,the rotational hardening limit surface mb, the initial degree ofanisotropy, the initial void ratio, the height of soilfoundation/specimen, the water surface level, the initial overburdenpressure on the surface of the soil foundation, the initialconsolidation pressure of the specimen, the cell pressure (tensile:positive), and settings of the unit system (time, force, and length).The input settings of soils are, however, not restricted to theseparameters. Only part of these parameters may be input as the settingsof soils or any other suitable parameters may be input as the settingsof soils. These settings of soils may be input in any suitable inputformat.

The soil-water coupled analyzer 20 of the embodiment inputs the settingsof the solid soil model as shown in the concrete example of FIG. 8. Theinput settings of the solid soil model in the file of FIG. 8 include thenumber of elements, the number of nodes, and the number of elements withapplication of distributed load rate in the solid soil model, the numberof nodes per element, the number of computation steps, and the number ofelastic elements on the 1^(st) line, each of node numbers and the typeof boundary condition in the x-direction, the value of the boundarycondition in the x-direction, the type of boundary condition in they-direction, and the value of the boundary condition in the y-directionwith regard to the corresponding node on the 2^(nd) to the 10^(th)lines, each of element numbers and the global node numbers mapped to thelocal node numbers with regard to the corresponding element on the11^(th) to the 14^(th) lines, each of node numbers and the coordinate inthe x-direction and the coordinate in the y-direction of thecorresponding node on the 15^(th) to the 23^(rd) lines, the number ofnodes requiring output of the reactive force and their node numbers onthe 24^(th) and 25^(th) lines. The input settings of the solid soilmodel in the file of FIG. 8 also include element numbers withapplication of distributed load rate and the side of application(internode of application), the value of the distributed load rate ofits left node in the x-direction, the value of the distributed load ratein its left node in the y-direction, the value of the distributed loadrate of its right node in the x-direction, and the value of thedistributed load rate of its right node in the y-direction with regardto the corresponding element on the 26^(th) and 27^(th) lines, and eachof element numbers and the element number adjacent to the internode 1 ofthe corresponding element, the element number adjacent to the internode2 of the corresponding element, the element number adjacent to theinternode 3 of the corresponding element, the element number adjacent tothe internode 4 of the corresponding element, and the soil type of thecorresponding element on the 28^(th) to the 31^(st) lines. The inputsettings of the solid soil model in the file of FIG. 8 further includethe number of conditions satisfying the fixed length between two nodeson the 32^(nd) line, the number of conditions satisfying the fixed anglebetween three nodes on the 33^(rd) line, the number of conditionssatisfying the fixed direction of the relative velocity of two nodesspecified by relative position vectors of former two nodes and lattertwo nodes on the 35^(th) line, the condition number satisfying the fixeddirection of the relative velocity of two nodes and specification ofnodes involved in the condition on the 36^(th) line, the number ofconditions satisfying the equal displacement (velocity) on the 37^(th)line, the number of regular waves input into the model and the number ofirregular waves input into the model on the 38^(th) line, the type ofregular wave vibration, the number of nodes under influence of theregular wave vibration, and the amplitude and the period of the regularwave vibration on the 39^(th) line, and specification of nodes underinfluence of the regular wave vibration on the 40^(th) line. The inputsettings of the solid soil model are, however, not restricted to theseparameters. Only part of these parameters may be input as the settingsof the solid soil model or any other suitable parameters may be input asthe settings of the soil soil model. These settings of the solid soilmodel may be input in any suitable input format.

The soil-water coupled analyzer 20 of the embodiment performs thecomputations with the pore water set to an uncompressible fluid. Thecomputations may alternatively be performed with the pore water set to acompressible fluid. In the latter case, the global tangent stiffnessequation (Equation (63)) is solved with some modification using a poreratio e estimated as the value of the element i at the time t=t+θΔt inthe course of convergence calculation and a density of pore water givenby Equation (93). In the coefficient matrix on the left side of Equation(63), (two) transposed L terms of the global L matrix (global volumechange rate matrix) are replaced by transposed sums of the correspondingL terms and a global Lc matrix (global pore water compressibilitymatrix) given by Equation (94).

$\begin{matrix}{\rho_{w} = {\rho_{w\; 0}{\exp\left( \frac{u - u_{w\; 0}}{K_{f}} \right)}}} & (93)\end{matrix}$ρ_(w0): density of pore water under reference water pressure u_(w0)K_(f): bulk modulus of elasticity of pore water

$\begin{matrix}{{L_{c} = {\begin{bmatrix}\begin{matrix}\begin{matrix}\begin{matrix}\begin{matrix}\begin{matrix}L_{c}^{e\; 1} \\L_{c}^{e\; 2}\end{matrix} \\\vdots\end{matrix} \\L_{c}^{e\; i}\end{matrix} \\\vdots\end{matrix} \\L_{c}^{e\;{({{NE} - 1})}}\end{matrix} \\L_{c}^{eNE}\end{bmatrix}\text{:}}}{{global}\mspace{14mu}{Lc}\mspace{14mu}{{matrix}{\mspace{14mu}\;}\left( {{global}\mspace{14mu}{pore}\mspace{14mu}{water}{compressibility}\mspace{14mu}{matrix}} \right)}}} & (94) \\{{L_{c}^{ei} = {- {\int_{{ve}^{i}}{\frac{\rho_{w}}{K_{f}}{\frac{e}{1 + e}\lbrack N\rbrack}^{T}\left( {{\lbrack N\rbrack\left\{ {\overset{.}{v}}^{N} \right\}} - \left\{ b \right\}} \right)\ {\mathbb{d}v}\text{:}}}}}{{element}\mspace{20mu}{Lc}\mspace{14mu}{{matrix}{\mspace{14mu}\;}\left( {{element}\mspace{14mu}{water}{compressibility}\mspace{14mu}{matrix}} \right)}\mspace{14mu}{of}\mspace{14mu}{element}\mspace{14mu} i}} & (95)\end{matrix}$

In the coefficient matrix on the left side of Equation (63), the globalH matrix (Equation (16)) is replaced by a modified global H matrix givenby Equation (96). Equation (98) is added to the lower term on the rightside of Equation (63).

$\begin{matrix}{{H_{c} = {\begin{bmatrix}\begin{matrix}\begin{matrix}\begin{matrix}\begin{matrix}\begin{matrix}H_{c}^{e\; 1} \\H_{c}^{e\; 2}\end{matrix} \\\vdots\end{matrix} \\H_{c}^{ei}\end{matrix} \\\vdots\end{matrix} \\H_{c}^{e\;{({{NE} - 1})}}\end{matrix} \\H_{c}^{eNE}\end{bmatrix}\text{:}}}{{modified}\mspace{14mu}{global}\mspace{14mu} H\mspace{14mu}{{matrix}{\mspace{14mu}\;}\left( {{modified}{global}\mspace{14mu}{water}\mspace{14mu}{permeability}\mspace{14mu}{matrix}} \right)}}} & (96) \\{{H_{c}^{ei} = {\left\lbrack {\alpha_{1}^{ei}\mspace{14mu}\ldots\mspace{14mu}\alpha_{m}^{ei}\mspace{14mu}\ldots\mspace{14mu}{\alpha_{s}^{ei}\left( {{- {\sum\limits_{m = 1}^{s}\alpha_{m}^{ei}}} - {\frac{2}{K_{f}}\frac{e}{1 + e}\frac{1}{{\theta\Delta}\; t}}} \right)}} \right\rbrack\text{:}}}{{element}\mspace{14mu}{Hc}\mspace{14mu}{{matrix}{\mspace{14mu}\;}\left( {{modified}\mspace{14mu}{element}\mspace{14mu}{water}\text{}{permeability}\mspace{14mu}{matrix}} \right)}\mspace{20mu}{of}\mspace{14mu}{element}\mspace{14mu} i}} & (97) \\{{\left\{ {{\overset{.}{f}}_{c}\left( {{\theta\Delta}\; t} \right)} \right\} ❘_{t}} = {\begin{Bmatrix}{{- 2}\frac{1}{K_{f}}\frac{e}{1 + e}\left( {u^{1} + {\frac{1}{2}{{\overset{.}{u}}^{1}\left( {{\theta\Delta}\; t} \right)}}} \right)} \\{{- 2}\frac{1}{K_{f}}\frac{e}{1 + e}\left( {u^{2} + {\frac{1}{2}{{\overset{.}{u}}^{2}\left( {{\theta\Delta}\; t} \right)}}} \right)} \\\vdots \\{{- 2}\frac{1}{K_{f}}\frac{e}{1 + e}\left( {u^{i} + {\frac{1}{2}{{\overset{.}{u}}^{i}\left( {{\theta\Delta}\; t} \right)}}} \right)} \\\vdots \\{{- 2}\frac{1}{K_{f}}\frac{e}{1 + e}\left( {u^{({{NE} - 1})} + {\frac{1}{2}{{\overset{.}{u}}^{({{NE} - 1})}\left( {{\theta\Delta}\; t} \right)}}} \right)} \\{{- 2}\frac{1}{K_{f}}\frac{e}{1 + e}\left( {u^{NE} + {\frac{1}{2}{{\overset{.}{u}}^{NE}\left( {{\theta\Delta}\; t} \right)}}} \right)}\end{Bmatrix}❘_{t}}} & (98)\end{matrix}$

In this modification, the soil-water coupled analysis program 30 adoptsa modified flow shown in the flowchart of FIG. 17, in place of theTamura/Christian model of the algorithm shown in the flowchart of FIG.2. The modified flow of FIG. 17 computes the modified global H matrix byEquation (96) (step S146) and the global pore water compressibilitymatrix (global Lc matrix) by Equation (94) (step S147), instead of thecomputation of the global H matrix at step S145 in the flow of FIG. 2.The modified flow also calculates Equation (98) in the computation ofthe total load increment vector and the total flow increment vector atstep S150. The density of pore water calculated from the water pressureat the element center according to Equation (93) is substituted into therelevant terms (for example, Equations (35) and (36)). The iterativeloop to detection of the convergence at step S190 is from step S140 tostep S180 in this modified flow. It is desirable to add a coefficient ofvolume compressibility Kf of pore water to the input settings of soilsshown as the concrete example of FIG. 7. The coefficient of permeabilitymay be varied corresponding to the void ratio.

FIG. 18 shows time-settlement curves after a start of earthquake motionimmediately below embankment as an example of computation withconsideration of the compressibility of pore water. The ‘earthquakemotion’ of about 100 gal at the maximum is given at the cycle ofapproximately 1 second from the bottom of a loose sand foundation(layered thickness of 20 m, embankment height of 4 m) having the varyingcompressibility of pore water only in a lower area of the embankmentloading. The time-settlement curve for the compressible pore water has asignificant compaction in the lower area of the embankment during theearthquake. The time-settlement curve for the compressible pore waterhas the greater depth of settlement in the initial stage after the endof the earthquake but eventually the less depth of settlement withelapse of time, compared with the time-settlement curve for theuncompressed pore water. The general conditions of the soil foundationare expected to change with the variation in compressibility of porewater in the lower area of the loading embankment. In this example ofcomputation, the density ρ_(w0) of compressible pore water at areference water pressure u_(w0) is varied to prevent the potentialchange of the conditions of the soil foundation.

The soil-water coupled analyzer 20 of the embodiment may execute thesoil-water coupled analysis program plural times with the results of theanalysis set for new analysis conditions or initial conditions. Theiterative execution of the soil-water coupled analysis program may belimited to a preset number of times, or execution of the soil-watercoupled analysis program may be allowed only once.

The soil-water coupled analyzer 20 of the embodiment computes theelement K matrix according to Equation (39) as the sum of the first termrepresenting the tangent stiffness of giving the mechanical behaviors ofthe specified soil materials (for example, sand and clay), the secondterm representing the contribution of the geometric time change of thesoil to the tangent stiffness, the third term representing thecontribution to the tangent stiffness under application of a significantacceleration to the soil skeleton, for example, under application of animpact load, and the fourth term representing the contribution to thetangent stiffness under application of a body force. According to theanalysis conditions, part of these terms may be omitted from thecomputation of the element K matrix.

The soil-water coupled analyzer 20 of the embodiment is applicable tothe analyses for any of the two-dimensional plane strain condition, theaxisymmetric condition, and the three-dimensional condition. Thesoil-water coupled analyzer may be designed to allow the analysis onlyfor the two-dimensional plane strain condition or the analyses only forthe two-dimensional plane strain condition and the axisymmetriccondition.

The soil-water coupled analyzer 20 of the embodiment is adopted toestimate the deformation behaviors (including bifurcation behavior andpost-peak behavior) and self weight consolidation analysis of thespecimen, the consolidation settlement and lateral movement behaviors(including ‘secondary consolidation’ long-term continuous settlementbehavior) of soil under loading like road embankment, the soildeformation behaviors in excavation, the liquefaction andliquefaction-induced settlement behaviors under repetitive loading likeearthquake, the interaction with an elastic structure (deformationbehaviors of the piled raft foundation), and the deformation behaviorsof reinforced soil. The soil-water coupled analyzer 20 of the embodimentis also adopted to study and investigate the water pumping-induced soilsubsidence, the soil bearing capacity, the stability/instability by thegroundwater flow, the compacted soil improvement with the cylindricalcavity expansion by sand piles, the soil improvement by improvement ofmass permeability, and the stability of an elasto-plastic embankmentstructure. In the case of such estimations and studies, execution of thesoil-water coupled analysis program plural times may be required withthe results of the analysis set for new analysis conditions or initialconditions. The soil-water coupled analyzer 20 may be adopted forestimation or investigation of only part or only one of these behaviorsand problems.

The embodiment regards the soil-water coupled analyzer 20 with thesoil-water coupled analysis program 30 installed therein. The techniqueof the invention may be attained by the soil-water coupled analysisprogram 30 or by the computation module 34 of the soil-water coupledanalysis program 30.

The embodiment and its modified examples discussed above are to beconsidered in all aspects as illustrative and not restrictive. There maybe many modifications, changes, and alterations without departing fromthe scope or spirit of the main characteristics of the presentinvention.

The present application claims priority from Japanese patent applicationNo. 2005-301415 filed on Oct. 17, 2005, the entire contents of which arehereby incorporated by reference into this application.

INDUSTRIAL APPLICABILITY

The technique of the invention is preferably applicable to themanufacturing industries utilizing the soil-water coupled analyzer.

DESCRIPTION OF NOTATIONS

-   l_(m): relative position vector from center of gravity of element i    to center of gravity of element m-   h: total water head at center of gravity of element i-   h_(m): total water head at center of gravity of element m-   h_(c): total water head at intersection of l_(i) and interface    between element i and element m-   ν′_(i) ^(m): flow rate vector of pore water in element i toward    element m-   ν′_(m) ^(m): flow rate vector of pore water in element m from    element i-   r_(a), r_(b): distance from axis of rotational symmetry to two nodes    at interface between element i and element m

1. A soil-water coupled analyzer comprising: a processor that: inputsdata on conditions of a soil skeleton and data on external forces;performs deformation analysis of the soil skeleton, wherein thedeformation analysis formulates simultaneous motion equations of ratetype with respect to a pore water pressure and a velocity of the solidskeleton including a jerk term of a second order derivative term on atime of a velocity vector of the soil skeleton, based on the input dataon the conditions of the soil skeleton and the input data on theexternal forces, and integrates the formulated simultaneous motionequations as needed under a given geometric boundary condition with theregard to any of the displacements, a displacement rate, and anacceleration, a given dynamic boundary condition with regard to eitherof a stress and a stress rate, and given hydraulic boundary conditionswith regard to either of the pore water pressure and a total water headand a flow of pore water, so as to obtain time history responses of avelocity field and pore water pressure field; and outputs the obtainedtime history responses.
 2. The soil-water coupled analyzer in accordancewith claim 1, wherein the analysis module calculates multiple termsbased on the input data on the conditions of the solid skeleton and theinput data on the external forces and formulates the simultaneous motionequations of rate type using the multiple calculated terms, where themultiple terms include at least one of the first term regarding a volumechange rate of the soil skeleton over time, a second term regarding awater permeability of soil, a third term regarding a mass, and a fourthterm regarding a tangent stiffness of the soil skeleton.
 3. Thesoil-water coupled analyzer in accordance with claim 2, wherein theanalysis module formulates the simultaneous motion equations of ratetype using all of the first term, the second term, the third term, andthe fourth term.
 4. The soil-water coupled analyzer in accordance withclaim 1, wherein the analysis module adopts a numerical analysistechnique, such as a finite element method or a difference method, toformulate the simultaneous motion equations of rate type.
 5. Thesoil-water coupled analyzer in accordance with claim 1, wherein the dataon the conditions of the soil skeleton include data on properties ofmultiple elements constituting the soil skeleton and data on relationsbetween adjacent elements.
 6. The soil-water coupled analyzer inaccordance with claim 5, wherein the data on the properties of themultiple elements include data on properties of soils.
 7. The soil-watercoupled analyzer in accordance with claim 1, wherein the data on theexternal forces include data on a load and a displacement applied to asoil-water coupled system.
 8. The soil-water coupled analyzer inaccordance with claim 1, wherein the data on the external forces includedata on a vibration applied to a soil-water coupled system.
 9. Thesoil-water coupled analyzer in accordance with claim 8, wherein theanalysis module performs the deformation analysis plural times withsetting a result of a current cycle of the analysis prior to the changeof the data on the external forces to an initial state of the data onthe conditions of the soil skeleton for a next cycle of the analysis.10. The soil-water coupled analyzer in accordance with claim 1, whereinthe analysis module performs the deformation analysis plural times withchanging the data on the external forces.
 11. A non-transitory computerreadable medium comprising programming for execution by a processor, theprogramming when executed, performing a soil-water coupled analysismethod, the method comprising: performing deformation analysis of a soilskeleton, wherein the deformation analysis: formulates simultaneousmotion equations of rate type with respect to a pore water pressure anda velocity of the soil skeleton including a jerk term of a second orderderivative term on a time of a velocity vector of the soil skeleton,based on data on conditions of the soil skeleton and data on externalforces, and integrates the formulated simultaneous motion equations asneeded under a given geometric boundary condition with regard to any ofa displacement, a displacement rate, and an acceleration, a givendynamic boundary condition with regard to either of a stress and astress rate, and given hydraulic boundary conditions with regard toeither of the pore water pressure and a total water head and a flow ofpore water, so as to obtain time history responses of a velocity fieldand a pore water pressure field.
 12. The computer readable medium inaccordance with claim 11, wherein the soil-water coupled analysis methodcalculates multiple terms based on the data on the conditions of thesoil skeleton and the data on the external forces and formulates thesimultaneous motion equations of rate type using the multiple calculatedterms, where the multiple terms include at least one of a first termregarding a volume change rate of the soil skeleton over time, a secondterm regarding a water permeability of soil, a third term regarding amass, and a fourth term regarding a tangent stiffness of the soilskeleton.